--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/array/dSparse.cc	Fri Aug 31 20:00:20 2012 -0700
@@ -0,0 +1,7997 @@
+/*
+
+Copyright (C) 2004-2012 David Bateman
+Copyright (C) 1998-2004 Andy Adler
+Copyright (C) 2010 VZLU Prague
+
+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 <cfloat>
+
+#include <iostream>
+#include <vector>
+#include <functional>
+
+#include "quit.h"
+#include "lo-ieee.h"
+#include "lo-mappers.h"
+#include "f77-fcn.h"
+#include "dRowVector.h"
+#include "oct-locbuf.h"
+
+#include "dDiagMatrix.h"
+#include "CSparse.h"
+#include "boolSparse.h"
+#include "dSparse.h"
+#include "functor.h"
+#include "oct-spparms.h"
+#include "SparsedbleLU.h"
+#include "MatrixType.h"
+#include "oct-sparse.h"
+#include "sparse-util.h"
+#include "SparsedbleCHOL.h"
+#include "SparseQR.h"
+
+#include "Sparse-diag-op-defs.h"
+
+#include "Sparse-perm-op-defs.h"
+
+// Define whether to use a basic QR solver or one that uses a Dulmange
+// Mendelsohn factorization to seperate the problem into under-determined,
+// well-determined and over-determined parts and solves them seperately
+#ifndef USE_QRSOLVE
+#include "sparse-dmsolve.cc"
+#endif
+
+// Fortran functions we call.
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (dgbtrf, DGBTRF) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*, const octave_idx_type&,
+                             octave_idx_type*, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (dgbtrs, DGBTRS) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const double*, const octave_idx_type&,
+                             const octave_idx_type*, double*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dgbcon, DGBCON) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, const octave_idx_type*,
+                             const double&, double&, double*,
+                             octave_idx_type*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dpbtrf, DPBTRF) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*, const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dpbtrs, DPBTRS) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dpbcon, DPBCON) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             double*, const octave_idx_type&,
+                             const double&, double&, double*,
+                             octave_idx_type*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+  F77_RET_T
+  F77_FUNC (dptsv, DPTSV) (const octave_idx_type&, const octave_idx_type&,
+                           double*, double*, double*, const octave_idx_type&,
+                           octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (dgtsv, DGTSV) (const octave_idx_type&, const octave_idx_type&,
+                           double*, double*, double*, double*,
+                           const octave_idx_type&, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (dgttrf, DGTTRF) (const octave_idx_type&, double*, double*,
+                             double*, double*, octave_idx_type*,
+                             octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (dgttrs, DGTTRS) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const double*, const double*, const double*,
+                             const double*, const octave_idx_type*,
+                             double *, const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zptsv, ZPTSV) (const octave_idx_type&, const octave_idx_type&,
+                           double*, Complex*, Complex*, const octave_idx_type&,
+                           octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (zgtsv, ZGTSV) (const octave_idx_type&, const octave_idx_type&,
+                           Complex*, Complex*, Complex*, Complex*,
+                           const octave_idx_type&, octave_idx_type&);
+
+}
+
+SparseMatrix::SparseMatrix (const SparseBoolMatrix &a)
+  : MSparse<double> (a.rows (), a.cols (), a.nnz ())
+{
+  octave_idx_type nc = cols ();
+  octave_idx_type nz = a.nnz ();
+
+  for (octave_idx_type i = 0; i < nc + 1; i++)
+    cidx (i) = a.cidx (i);
+
+  for (octave_idx_type i = 0; i < nz; i++)
+    {
+      data (i) = a.data (i);
+      ridx (i) = a.ridx (i);
+    }
+}
+
+SparseMatrix::SparseMatrix (const DiagMatrix& a)
+  : MSparse<double> (a.rows (), a.cols (), a.length ())
+{
+  octave_idx_type j = 0, l = a.length ();
+  for (octave_idx_type i = 0; i < l; i++)
+    {
+      cidx (i) = j;
+      if (a(i, i) != 0.0)
+        {
+          data (j) = a(i, i);
+          ridx (j) = i;
+          j++;
+        }
+    }
+  for (octave_idx_type i = l; i <= a.cols (); i++)
+    cidx (i) = j;
+}
+
+bool
+SparseMatrix::operator == (const SparseMatrix& a) const
+{
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nz = nnz ();
+  octave_idx_type nr_a = a.rows ();
+  octave_idx_type nc_a = a.cols ();
+  octave_idx_type nz_a = a.nnz ();
+
+  if (nr != nr_a || nc != nc_a || nz != nz_a)
+    return false;
+
+  for (octave_idx_type i = 0; i < nc + 1; i++)
+    if (cidx (i) != a.cidx (i))
+        return false;
+
+  for (octave_idx_type i = 0; i < nz; i++)
+    if (data (i) != a.data (i) || ridx (i) != a.ridx (i))
+      return false;
+
+  return true;
+}
+
+bool
+SparseMatrix::operator != (const SparseMatrix& a) const
+{
+  return !(*this == a);
+}
+
+bool
+SparseMatrix::is_symmetric (void) const
+{
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+
+  if (nr == nc && nr > 0)
+    {
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+            {
+              octave_idx_type ri = ridx (i);
+
+              if (ri != j)
+                {
+                  bool found = false;
+
+                  for (octave_idx_type k = cidx (ri); k < cidx (ri+1); k++)
+                    {
+                      if (ridx (k) == j)
+                        {
+                          if (data (i) == data (k))
+                            found = true;
+                          break;
+                        }
+                    }
+
+                  if (! found)
+                    return false;
+                }
+            }
+        }
+
+      return true;
+    }
+
+  return false;
+}
+
+SparseMatrix&
+SparseMatrix::insert (const SparseMatrix& a, octave_idx_type r, octave_idx_type c)
+{
+  MSparse<double>::insert (a, r, c);
+  return *this;
+}
+
+SparseMatrix&
+SparseMatrix::insert (const SparseMatrix& a, const Array<octave_idx_type>& indx)
+{
+  MSparse<double>::insert (a, indx);
+  return *this;
+}
+
+SparseMatrix
+SparseMatrix::max (int dim) const
+{
+  Array<octave_idx_type> dummy_idx;
+  return max (dummy_idx, dim);
+}
+
+SparseMatrix
+SparseMatrix::max (Array<octave_idx_type>& idx_arg, int dim) const
+{
+  SparseMatrix result;
+  dim_vector dv = dims ();
+
+  if (dv.numel () == 0 || dim >= dv.length ())
+    return result;
+
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
+  octave_idx_type nr = dv(0);
+  octave_idx_type nc = dv(1);
+
+  if (dim == 0)
+    {
+      idx_arg.clear (1, nc);
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          double tmp_max = octave_NaN;
+          octave_idx_type idx_j = 0;
+          for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+            {
+              if (ridx (i) != idx_j)
+                break;
+              else
+                idx_j++;
+            }
+
+          if (idx_j != nr)
+            tmp_max = 0.;
+
+          for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+            {
+              double tmp = data (i);
+
+              if (xisnan (tmp))
+                continue;
+              else if (xisnan (tmp_max) || tmp > tmp_max)
+                {
+                  idx_j = ridx (i);
+                  tmp_max = tmp;
+                }
+
+            }
+
+          idx_arg.elem (j) = xisnan (tmp_max) ? 0 : idx_j;
+          if (tmp_max != 0.)
+            nel++;
+        }
+
+      result = SparseMatrix (1, nc, nel);
+
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          double tmp = elem (idx_arg(j), j);
+          if (tmp != 0.)
+            {
+              result.xdata (ii) = tmp;
+              result.xridx (ii++) = 0;
+            }
+          result.xcidx (j+1) = ii;
+
+        }
+    }
+  else
+    {
+      idx_arg.resize (dim_vector  (nr, 1), 0);
+
+      for (octave_idx_type i = cidx (0); i < cidx (1); i++)
+        idx_arg.elem (ridx (i)) = -1;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            if (idx_arg.elem (i) != -1)
+              continue;
+            bool found = false;
+            for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
+              if (ridx (k) == i)
+                {
+                  found = true;
+                  break;
+                }
+
+            if (!found)
+              idx_arg.elem (i) = j;
+
+          }
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+            {
+              octave_idx_type ir = ridx (i);
+              octave_idx_type ix = idx_arg.elem (ir);
+              double tmp = data (i);
+
+              if (xisnan (tmp))
+                continue;
+              else if (ix == -1 || tmp > elem (ir, ix))
+                idx_arg.elem (ir) = j;
+            }
+        }
+
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nr; j++)
+        if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
+          nel++;
+
+      result = SparseMatrix (nr, 1, nel);
+
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+      result.xcidx (1) = nel;
+      for (octave_idx_type j = 0; j < nr; j++)
+        {
+          if (idx_arg(j) == -1)
+            {
+              idx_arg(j) = 0;
+              result.xdata (ii) = octave_NaN;
+              result.xridx (ii++) = j;
+            }
+          else
+            {
+              double tmp = elem (j, idx_arg(j));
+              if (tmp != 0.)
+                {
+                  result.xdata (ii) = tmp;
+                  result.xridx (ii++) = j;
+                }
+            }
+        }
+    }
+
+  return result;
+}
+
+SparseMatrix
+SparseMatrix::min (int dim) const
+{
+  Array<octave_idx_type> dummy_idx;
+  return min (dummy_idx, dim);
+}
+
+SparseMatrix
+SparseMatrix::min (Array<octave_idx_type>& idx_arg, int dim) const
+{
+  SparseMatrix result;
+  dim_vector dv = dims ();
+
+  if (dv.numel () == 0 || dim >= dv.length ())
+    return result;
+
+  if (dim < 0)
+    dim = dv.first_non_singleton ();
+
+  octave_idx_type nr = dv(0);
+  octave_idx_type nc = dv(1);
+
+  if (dim == 0)
+    {
+      idx_arg.clear (1, nc);
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          double tmp_min = octave_NaN;
+          octave_idx_type idx_j = 0;
+          for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+            {
+              if (ridx (i) != idx_j)
+                break;
+              else
+                idx_j++;
+            }
+
+          if (idx_j != nr)
+            tmp_min = 0.;
+
+          for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+            {
+              double tmp = data (i);
+
+              if (xisnan (tmp))
+                continue;
+              else if (xisnan (tmp_min) || tmp < tmp_min)
+                {
+                  idx_j = ridx (i);
+                  tmp_min = tmp;
+                }
+
+            }
+
+          idx_arg.elem (j) = xisnan (tmp_min) ? 0 : idx_j;
+          if (tmp_min != 0.)
+            nel++;
+        }
+
+      result = SparseMatrix (1, nc, nel);
+
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          double tmp = elem (idx_arg(j), j);
+          if (tmp != 0.)
+            {
+              result.xdata (ii) = tmp;
+              result.xridx (ii++) = 0;
+            }
+          result.xcidx (j+1) = ii;
+
+        }
+    }
+  else
+    {
+      idx_arg.resize (dim_vector (nr, 1), 0);
+
+      for (octave_idx_type i = cidx (0); i < cidx (1); i++)
+        idx_arg.elem (ridx (i)) = -1;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = 0; i < nr; i++)
+          {
+            if (idx_arg.elem (i) != -1)
+              continue;
+            bool found = false;
+            for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
+              if (ridx (k) == i)
+                {
+                  found = true;
+                  break;
+                }
+
+            if (!found)
+              idx_arg.elem (i) = j;
+
+          }
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+            {
+              octave_idx_type ir = ridx (i);
+              octave_idx_type ix = idx_arg.elem (ir);
+              double tmp = data (i);
+
+              if (xisnan (tmp))
+                continue;
+              else if (ix == -1 || tmp < elem (ir, ix))
+                idx_arg.elem (ir) = j;
+            }
+        }
+
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nr; j++)
+        if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
+          nel++;
+
+      result = SparseMatrix (nr, 1, nel);
+
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+      result.xcidx (1) = nel;
+      for (octave_idx_type j = 0; j < nr; j++)
+        {
+          if (idx_arg(j) == -1)
+            {
+              idx_arg(j) = 0;
+              result.xdata (ii) = octave_NaN;
+              result.xridx (ii++) = j;
+            }
+          else
+            {
+              double tmp = elem (j, idx_arg(j));
+              if (tmp != 0.)
+                {
+                  result.xdata (ii) = tmp;
+                  result.xridx (ii++) = j;
+                }
+            }
+        }
+    }
+
+  return result;
+}
+
+RowVector
+SparseMatrix::row (octave_idx_type i) const
+{
+  octave_idx_type nc = columns ();
+  RowVector retval (nc, 0);
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
+      {
+        if (ridx (k) == i)
+          {
+            retval(j) = data (k);
+            break;
+          }
+      }
+
+  return retval;
+}
+
+ColumnVector
+SparseMatrix::column (octave_idx_type i) const
+{
+  octave_idx_type nr = rows ();
+  ColumnVector retval (nr, 0);
+
+  for (octave_idx_type k = cidx (i); k < cidx (i+1); k++)
+    retval(ridx (k)) = data (k);
+
+  return retval;
+}
+
+SparseMatrix
+SparseMatrix::concat (const SparseMatrix& rb, const Array<octave_idx_type>& ra_idx)
+{
+  // Don't use numel to avoid all possiblity of an overflow
+  if (rb.rows () > 0 && rb.cols () > 0)
+    insert (rb, ra_idx(0), ra_idx(1));
+  return *this;
+}
+
+SparseComplexMatrix
+SparseMatrix::concat (const SparseComplexMatrix& rb, const Array<octave_idx_type>& ra_idx)
+{
+  SparseComplexMatrix retval (*this);
+  if (rb.rows () > 0 && rb.cols () > 0)
+    retval.insert (rb, ra_idx(0), ra_idx(1));
+  return retval;
+}
+
+SparseMatrix
+real (const SparseComplexMatrix& a)
+{
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+  octave_idx_type nz = a.nnz ();
+  SparseMatrix r (nr, nc, nz);
+
+  for (octave_idx_type i = 0; i < nc +1; i++)
+    r.cidx (i) = a.cidx (i);
+
+  for (octave_idx_type i = 0; i < nz; i++)
+    {
+      r.data (i) = std::real (a.data (i));
+      r.ridx (i) = a.ridx (i);
+    }
+
+  return r;
+}
+
+SparseMatrix
+imag (const SparseComplexMatrix& a)
+{
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+  octave_idx_type nz = a.nnz ();
+  SparseMatrix r (nr, nc, nz);
+
+  for (octave_idx_type i = 0; i < nc +1; i++)
+    r.cidx (i) = a.cidx (i);
+
+  for (octave_idx_type i = 0; i < nz; i++)
+    {
+      r.data (i) = std::imag (a.data (i));
+      r.ridx (i) = a.ridx (i);
+    }
+
+  return r;
+}
+
+SparseMatrix
+atan2 (const double& x, const SparseMatrix& y)
+{
+  octave_idx_type nr = y.rows ();
+  octave_idx_type nc = y.cols ();
+
+  if (x == 0.)
+    return SparseMatrix (nr, nc);
+  else
+    {
+      // Its going to be basically full, so this is probably the
+      // best way to handle it.
+      Matrix tmp (nr, nc, atan2 (x, 0.));
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = y.cidx (j); i < y.cidx (j+1); i++)
+          tmp.elem (y.ridx (i), j) = atan2 (x, y.data (i));
+
+      return SparseMatrix (tmp);
+    }
+}
+
+SparseMatrix
+atan2 (const SparseMatrix& x, const double& y)
+{
+  octave_idx_type nr = x.rows ();
+  octave_idx_type nc = x.cols ();
+  octave_idx_type nz = x.nnz ();
+
+  SparseMatrix retval (nr, nc, nz);
+
+  octave_idx_type ii = 0;
+  retval.xcidx (0) = 0;
+  for (octave_idx_type i = 0; i < nc; i++)
+    {
+      for (octave_idx_type j = x.cidx (i); j < x.cidx (i+1); j++)
+        {
+          double tmp = atan2 (x.data (j), y);
+          if (tmp != 0.)
+            {
+              retval.xdata (ii) = tmp;
+              retval.xridx (ii++) = x.ridx (j);
+            }
+        }
+      retval.xcidx (i+1) = ii;
+    }
+
+  if (ii != nz)
+    {
+      SparseMatrix retval2 (nr, nc, ii);
+      for (octave_idx_type i = 0; i < nc+1; i++)
+        retval2.xcidx (i) = retval.cidx (i);
+      for (octave_idx_type i = 0; i < ii; i++)
+        {
+          retval2.xdata (i) = retval.data (i);
+          retval2.xridx (i) = retval.ridx (i);
+        }
+      return retval2;
+    }
+  else
+    return retval;
+}
+
+SparseMatrix
+atan2 (const SparseMatrix& x, const SparseMatrix& y)
+{
+  SparseMatrix r;
+
+  if ((x.rows () == y.rows ()) && (x.cols () == y.cols ()))
+    {
+      octave_idx_type x_nr = x.rows ();
+      octave_idx_type x_nc = x.cols ();
+
+      octave_idx_type y_nr = y.rows ();
+      octave_idx_type y_nc = y.cols ();
+
+      if (x_nr != y_nr || x_nc != y_nc)
+        gripe_nonconformant ("atan2", x_nr, x_nc, y_nr, y_nc);
+      else
+        {
+          r = SparseMatrix (x_nr, x_nc, (x.nnz () + y.nnz ()));
+
+          octave_idx_type jx = 0;
+          r.cidx (0) = 0;
+          for (octave_idx_type i = 0 ; i < x_nc ; i++)
+            {
+              octave_idx_type  ja = x.cidx (i);
+              octave_idx_type  ja_max = x.cidx (i+1);
+              bool ja_lt_max= ja < ja_max;
+
+              octave_idx_type  jb = y.cidx (i);
+              octave_idx_type  jb_max = y.cidx (i+1);
+              bool jb_lt_max = jb < jb_max;
+
+              while (ja_lt_max || jb_lt_max )
+                {
+                  octave_quit ();
+                  if ((! jb_lt_max) ||
+                      (ja_lt_max && (x.ridx (ja) < y.ridx (jb))))
+                    {
+                      r.ridx (jx) = x.ridx (ja);
+                      r.data (jx) = atan2 (x.data (ja), 0.);
+                      jx++;
+                      ja++;
+                      ja_lt_max= ja < ja_max;
+                    }
+                  else if (( !ja_lt_max ) ||
+                           (jb_lt_max && (y.ridx (jb) < x.ridx (ja)) ) )
+                    {
+                      jb++;
+                      jb_lt_max= jb < jb_max;
+                    }
+                  else
+                    {
+                      double tmp = atan2 (x.data (ja), y.data (jb));
+                      if (tmp != 0.)
+                        {
+                          r.data (jx) = tmp;
+                          r.ridx (jx) = x.ridx (ja);
+                          jx++;
+                        }
+                      ja++;
+                      ja_lt_max= ja < ja_max;
+                      jb++;
+                      jb_lt_max= jb < jb_max;
+                    }
+                }
+              r.cidx (i+1) = jx;
+            }
+
+          r.maybe_compress ();
+        }
+    }
+  else
+    (*current_liboctave_error_handler) ("matrix size mismatch");
+
+  return r;
+}
+
+SparseMatrix
+SparseMatrix::inverse (void) const
+{
+  octave_idx_type info;
+  double rcond;
+  MatrixType mattype (*this);
+  return inverse (mattype, info, rcond, 0, 0);
+}
+
+SparseMatrix
+SparseMatrix::inverse (MatrixType& mattype) const
+{
+  octave_idx_type info;
+  double rcond;
+  return inverse (mattype, info, rcond, 0, 0);
+}
+
+SparseMatrix
+SparseMatrix::inverse (MatrixType& mattype, octave_idx_type& info) const
+{
+  double rcond;
+  return inverse (mattype, info, rcond, 0, 0);
+}
+
+SparseMatrix
+SparseMatrix::dinverse (MatrixType &mattyp, octave_idx_type& info,
+                        double& rcond, const bool,
+                        const bool calccond) const
+{
+  SparseMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  info = 0;
+
+  if (nr == 0 || nc == 0 || nr != nc)
+    (*current_liboctave_error_handler) ("inverse requires square matrix");
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattyp.type ();
+      mattyp.info ();
+
+      if (typ == MatrixType::Diagonal ||
+          typ == MatrixType::Permuted_Diagonal)
+        {
+          if (typ == MatrixType::Permuted_Diagonal)
+            retval = transpose ();
+          else
+            retval = *this;
+
+          // Force make_unique to be called
+          double *v = retval.data ();
+
+          if (calccond)
+            {
+              double dmax = 0., dmin = octave_Inf;
+              for (octave_idx_type i = 0; i < nr; i++)
+                {
+                  double tmp = fabs (v[i]);
+                  if (tmp > dmax)
+                    dmax = tmp;
+                  if (tmp < dmin)
+                    dmin = tmp;
+                }
+              rcond = dmin / dmax;
+            }
+
+          for (octave_idx_type i = 0; i < nr; i++)
+            v[i] = 1.0 / v[i];
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseMatrix
+SparseMatrix::tinverse (MatrixType &mattyp, octave_idx_type& info,
+                        double& rcond, const bool,
+                        const bool calccond) const
+{
+  SparseMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  info = 0;
+
+  if (nr == 0 || nc == 0 || nr != nc)
+    (*current_liboctave_error_handler) ("inverse requires square matrix");
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattyp.type ();
+      mattyp.info ();
+
+      if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper ||
+          typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
+        {
+          double anorm = 0.;
+          double ainvnorm = 0.;
+
+          if (calccond)
+            {
+              // Calculate the 1-norm of matrix for rcond calculation
+              for (octave_idx_type j = 0; j < nr; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          if (typ == MatrixType::Upper || typ == MatrixType::Lower)
+            {
+              octave_idx_type nz = nnz ();
+              octave_idx_type cx = 0;
+              octave_idx_type nz2 = nz;
+              retval = SparseMatrix (nr, nc, nz2);
+
+              for (octave_idx_type i = 0; i < nr; i++)
+                {
+                  octave_quit ();
+                  // place the 1 in the identity position
+                  octave_idx_type cx_colstart = cx;
+
+                  if (cx == nz2)
+                    {
+                      nz2 *= 2;
+                      retval.change_capacity (nz2);
+                    }
+
+                  retval.xcidx (i) = cx;
+                  retval.xridx (cx) = i;
+                  retval.xdata (cx) = 1.0;
+                  cx++;
+
+                  // iterate accross columns of input matrix
+                  for (octave_idx_type j = i+1; j < nr; j++)
+                    {
+                      double v = 0.;
+                      // iterate to calculate sum
+                      octave_idx_type colXp = retval.xcidx (i);
+                      octave_idx_type colUp = cidx (j);
+                      octave_idx_type rpX, rpU;
+
+                      if (cidx (j) == cidx (j+1))
+                        {
+                          (*current_liboctave_error_handler)
+                            ("division by zero");
+                          goto inverse_singular;
+                        }
+
+                      do
+                        {
+                          octave_quit ();
+                          rpX = retval.xridx (colXp);
+                          rpU = ridx (colUp);
+
+                          if (rpX < rpU)
+                            colXp++;
+                          else if (rpX > rpU)
+                            colUp++;
+                          else
+                            {
+                              v -= retval.xdata (colXp) * data (colUp);
+                              colXp++;
+                              colUp++;
+                            }
+                        } while ((rpX<j) && (rpU<j) &&
+                                 (colXp<cx) && (colUp<nz));
+
+                      // get A(m,m)
+                      if (typ == MatrixType::Upper)
+                        colUp = cidx (j+1) - 1;
+                      else
+                        colUp = cidx (j);
+                      double pivot = data (colUp);
+                      if (pivot == 0. || ridx (colUp) != j)
+                        {
+                          (*current_liboctave_error_handler)
+                            ("division by zero");
+                          goto inverse_singular;
+                        }
+
+                      if (v != 0.)
+                        {
+                          if (cx == nz2)
+                            {
+                              nz2 *= 2;
+                              retval.change_capacity (nz2);
+                            }
+
+                          retval.xridx (cx) = j;
+                          retval.xdata (cx) = v / pivot;
+                          cx++;
+                        }
+                    }
+
+                  // get A(m,m)
+                  octave_idx_type colUp;
+                  if (typ == MatrixType::Upper)
+                    colUp = cidx (i+1) - 1;
+                  else
+                    colUp = cidx (i);
+                  double pivot = data (colUp);
+                  if (pivot == 0. || ridx (colUp) != i)
+                    {
+                      (*current_liboctave_error_handler) ("division by zero");
+                      goto inverse_singular;
+                    }
+
+                  if (pivot != 1.0)
+                    for (octave_idx_type j = cx_colstart; j < cx; j++)
+                      retval.xdata (j) /= pivot;
+                }
+              retval.xcidx (nr) = cx;
+              retval.maybe_compress ();
+            }
+          else
+            {
+              octave_idx_type nz = nnz ();
+              octave_idx_type cx = 0;
+              octave_idx_type nz2 = nz;
+              retval = SparseMatrix (nr, nc, nz2);
+
+              OCTAVE_LOCAL_BUFFER (double, work, nr);
+              OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr);
+
+              octave_idx_type *perm = mattyp.triangular_perm ();
+              if (typ == MatrixType::Permuted_Upper)
+                {
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    rperm[perm[i]] = i;
+                }
+              else
+                {
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    rperm[i] = perm[i];
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    perm[rperm[i]] = i;
+                }
+
+              for (octave_idx_type i = 0; i < nr; i++)
+                {
+                  octave_quit ();
+                  octave_idx_type iidx = rperm[i];
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    work[j] = 0.;
+
+                  // place the 1 in the identity position
+                  work[iidx] = 1.0;
+
+                  // iterate accross columns of input matrix
+                  for (octave_idx_type j = iidx+1; j < nr; j++)
+                    {
+                      double v = 0.;
+                      octave_idx_type jidx = perm[j];
+                      // iterate to calculate sum
+                      for (octave_idx_type k = cidx (jidx);
+                           k < cidx (jidx+1); k++)
+                        {
+                          octave_quit ();
+                          v -= work[ridx (k)] * data (k);
+                        }
+
+                      // get A(m,m)
+                      double pivot;
+                      if (typ == MatrixType::Permuted_Upper)
+                        pivot = data (cidx (jidx+1) - 1);
+                      else
+                        pivot = data (cidx (jidx));
+                      if (pivot == 0.)
+                        {
+                          (*current_liboctave_error_handler)
+                            ("division by zero");
+                          goto inverse_singular;
+                        }
+
+                      work[j] = v / pivot;
+                    }
+
+                  // get A(m,m)
+                  octave_idx_type colUp;
+                  if (typ == MatrixType::Permuted_Upper)
+                    colUp = cidx (perm[iidx]+1) - 1;
+                  else
+                    colUp = cidx (perm[iidx]);
+
+                  double pivot = data (colUp);
+                  if (pivot == 0.)
+                    {
+                      (*current_liboctave_error_handler)
+                        ("division by zero");
+                      goto inverse_singular;
+                    }
+
+                  octave_idx_type new_cx = cx;
+                  for (octave_idx_type j = iidx; j < nr; j++)
+                    if (work[j] != 0.0)
+                      {
+                        new_cx++;
+                        if (pivot != 1.0)
+                          work[j] /= pivot;
+                      }
+
+                  if (cx < new_cx)
+                    {
+                      nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
+                      retval.change_capacity (nz2);
+                    }
+
+                  retval.xcidx (i) = cx;
+                  for (octave_idx_type j = iidx; j < nr; j++)
+                    if (work[j] != 0.)
+                      {
+                        retval.xridx (cx) = j;
+                        retval.xdata (cx++) = work[j];
+                      }
+                }
+
+              retval.xcidx (nr) = cx;
+              retval.maybe_compress ();
+            }
+
+          if (calccond)
+            {
+              // Calculate the 1-norm of inverse matrix for rcond calculation
+              for (octave_idx_type j = 0; j < nr; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = retval.cidx (j);
+                       i < retval.cidx (j+1); i++)
+                    atmp += fabs (retval.data (i));
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
+                }
+
+              rcond = 1. / ainvnorm / anorm;
+            }
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+
+ inverse_singular:
+  return SparseMatrix ();
+}
+
+SparseMatrix
+SparseMatrix::inverse (MatrixType &mattype, octave_idx_type& info,
+                       double& rcond, int, int calc_cond) const
+{
+  int typ = mattype.type (false);
+  SparseMatrix ret;
+
+  if (typ == MatrixType::Unknown)
+    typ = mattype.type (*this);
+
+  if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
+    ret = dinverse (mattype, info, rcond, true, calc_cond);
+  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
+    ret = tinverse (mattype, info, rcond, true, calc_cond).transpose ();
+  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
+    {
+      MatrixType newtype = mattype.transpose ();
+      ret = transpose ().tinverse (newtype, info, rcond, true, calc_cond);
+    }
+  else
+    {
+      if (mattype.is_hermitian ())
+        {
+          MatrixType tmp_typ (MatrixType::Upper);
+          SparseCHOL fact (*this, info, false);
+          rcond = fact.rcond ();
+          if (info == 0)
+            {
+              double rcond2;
+              SparseMatrix Q = fact.Q ();
+              SparseMatrix InvL = fact.L ().transpose ().tinverse (tmp_typ,
+                                           info, rcond2, true, false);
+              ret = Q * InvL.transpose () * InvL * Q.transpose ();
+            }
+          else
+            {
+              // Matrix is either singular or not positive definite
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Full;
+            }
+        }
+
+      if (!mattype.is_hermitian ())
+        {
+          octave_idx_type n = rows ();
+          ColumnVector Qinit(n);
+          for (octave_idx_type i = 0; i < n; i++)
+            Qinit(i) = i;
+
+          MatrixType tmp_typ (MatrixType::Upper);
+          SparseLU fact (*this, Qinit, Matrix (), false, false);
+          rcond = fact.rcond ();
+          double rcond2;
+          SparseMatrix InvL = fact.L ().transpose ().tinverse (tmp_typ,
+                                           info, rcond2, true, false);
+          SparseMatrix InvU = fact.U ().tinverse (tmp_typ, info, rcond2,
+                                           true, false).transpose ();
+          ret = fact.Pc ().transpose () * InvU * InvL * fact.Pr ();
+        }
+    }
+
+  return ret;
+}
+
+DET
+SparseMatrix::determinant (void) const
+{
+  octave_idx_type info;
+  double rcond;
+  return determinant (info, rcond, 0);
+}
+
+DET
+SparseMatrix::determinant (octave_idx_type& info) const
+{
+  double rcond;
+  return determinant (info, rcond, 0);
+}
+
+DET
+SparseMatrix::determinant (octave_idx_type& err, double& rcond, int) const
+{
+  DET retval;
+
+#ifdef HAVE_UMFPACK
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+
+  if (nr == 0 || nc == 0 || nr != nc)
+    {
+      retval = DET (1.0);
+    }
+  else
+    {
+      err = 0;
+
+      // Setup the control parameters
+      Matrix Control (UMFPACK_CONTROL, 1);
+      double *control = Control.fortran_vec ();
+      UMFPACK_DNAME (defaults) (control);
+
+      double tmp = octave_sparse_params::get_key ("spumoni");
+      if (!xisnan (tmp))
+        Control (UMFPACK_PRL) = tmp;
+
+      tmp = octave_sparse_params::get_key ("piv_tol");
+      if (!xisnan (tmp))
+        {
+          Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
+          Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
+        }
+
+      // Set whether we are allowed to modify Q or not
+      tmp = octave_sparse_params::get_key ("autoamd");
+      if (!xisnan (tmp))
+        Control (UMFPACK_FIXQ) = tmp;
+
+      // Turn-off UMFPACK scaling for LU
+      Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
+
+      UMFPACK_DNAME (report_control) (control);
+
+      const octave_idx_type *Ap = cidx ();
+      const octave_idx_type *Ai = ridx ();
+      const double *Ax = data ();
+
+      UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1, control);
+
+      void *Symbolic;
+      Matrix Info (1, UMFPACK_INFO);
+      double *info = Info.fortran_vec ();
+      int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai,
+                                         Ax, 0, &Symbolic, control, info);
+
+      if (status < 0)
+        {
+          (*current_liboctave_error_handler)
+            ("SparseMatrix::determinant symbolic factorization failed");
+
+          UMFPACK_DNAME (report_status) (control, status);
+          UMFPACK_DNAME (report_info) (control, info);
+
+          UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
+        }
+      else
+        {
+          UMFPACK_DNAME (report_symbolic) (Symbolic, control);
+
+          void *Numeric;
+          status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic,
+                                       &Numeric, control, info) ;
+          UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
+
+          rcond = Info (UMFPACK_RCOND);
+
+          if (status < 0)
+            {
+              (*current_liboctave_error_handler)
+                ("SparseMatrix::determinant numeric factorization failed");
+
+              UMFPACK_DNAME (report_status) (control, status);
+              UMFPACK_DNAME (report_info) (control, info);
+
+              UMFPACK_DNAME (free_numeric) (&Numeric);
+            }
+          else
+            {
+              UMFPACK_DNAME (report_numeric) (Numeric, control);
+
+              double c10, e10;
+
+              status = UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric, info);
+
+              if (status < 0)
+                {
+                  (*current_liboctave_error_handler)
+                    ("SparseMatrix::determinant error calculating determinant");
+
+                  UMFPACK_DNAME (report_status) (control, status);
+                  UMFPACK_DNAME (report_info) (control, info);
+                }
+              else
+                retval = DET (c10, e10, 10);
+
+              UMFPACK_DNAME (free_numeric) (&Numeric);
+            }
+        }
+    }
+#else
+  (*current_liboctave_error_handler) ("UMFPACK not installed");
+#endif
+
+  return retval;
+}
+
+Matrix
+SparseMatrix::dsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& err,
+                      double& rcond, solve_singularity_handler,
+                      bool calc_cond) const
+{
+  Matrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc < nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = Matrix (nc, b.cols (), 0.0);
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Diagonal ||
+          typ == MatrixType::Permuted_Diagonal)
+        {
+          retval.resize (nc, b.cols (), 0.);
+          if (typ == MatrixType::Diagonal)
+            for (octave_idx_type j = 0; j < b.cols (); j++)
+              for (octave_idx_type i = 0; i < nm; i++)
+                retval(i,j) = b(i,j) / data (i);
+          else
+            for (octave_idx_type j = 0; j < b.cols (); j++)
+              for (octave_idx_type k = 0; k < nc; k++)
+                for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                  retval(k,j) = b(ridx (i),j) / data (i);
+
+          if (calc_cond)
+            {
+              double dmax = 0., dmin = octave_Inf;
+              for (octave_idx_type i = 0; i < nm; i++)
+                {
+                  double tmp = fabs (data (i));
+                  if (tmp > dmax)
+                    dmax = tmp;
+                  if (tmp < dmin)
+                    dmin = tmp;
+                }
+              rcond = dmin / dmax;
+            }
+          else
+            rcond = 1.;
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseMatrix
+SparseMatrix::dsolve (MatrixType &mattype, const SparseMatrix& b,
+                      octave_idx_type& err, double& rcond,
+                      solve_singularity_handler, bool calc_cond) const
+{
+  SparseMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc < nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = SparseMatrix (nc, b.cols ());
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Diagonal ||
+          typ == MatrixType::Permuted_Diagonal)
+        {
+          octave_idx_type b_nc = b.cols ();
+          octave_idx_type b_nz = b.nnz ();
+          retval = SparseMatrix (nc, b_nc, b_nz);
+
+          retval.xcidx (0) = 0;
+          octave_idx_type ii = 0;
+          if (typ == MatrixType::Diagonal)
+            for (octave_idx_type j = 0; j < b_nc; j++)
+              {
+                for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                  {
+                    if (b.ridx (i) >= nm)
+                      break;
+                    retval.xridx (ii) = b.ridx (i);
+                    retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
+                  }
+                retval.xcidx (j+1) = ii;
+              }
+          else
+            for (octave_idx_type j = 0; j < b_nc; j++)
+              {
+                for (octave_idx_type l = 0; l < nc; l++)
+                  for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
+                    {
+                      bool found = false;
+                      octave_idx_type k;
+                      for (k = b.cidx (j); k < b.cidx (j+1); k++)
+                        if (ridx (i) == b.ridx (k))
+                          {
+                            found = true;
+                            break;
+                          }
+                      if (found)
+                        {
+                          retval.xridx (ii) = l;
+                          retval.xdata (ii++) = b.data (k) / data (i);
+                        }
+                    }
+                retval.xcidx (j+1) = ii;
+              }
+
+          if (calc_cond)
+            {
+              double dmax = 0., dmin = octave_Inf;
+              for (octave_idx_type i = 0; i < nm; i++)
+                {
+                  double tmp = fabs (data (i));
+                  if (tmp > dmax)
+                    dmax = tmp;
+                  if (tmp < dmin)
+                    dmin = tmp;
+                }
+              rcond = dmin / dmax;
+            }
+          else
+            rcond = 1.;
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+ComplexMatrix
+SparseMatrix::dsolve (MatrixType &mattype, const ComplexMatrix& b,
+                      octave_idx_type& err, double& rcond,
+                      solve_singularity_handler, bool calc_cond) const
+{
+  ComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc < nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Diagonal ||
+          typ == MatrixType::Permuted_Diagonal)
+        {
+          retval.resize (nc, b.cols (), 0);
+          if (typ == MatrixType::Diagonal)
+            for (octave_idx_type j = 0; j < b.cols (); j++)
+                for (octave_idx_type i = 0; i < nm; i++)
+                  retval(i,j) = b(i,j) / data (i);
+          else
+            for (octave_idx_type j = 0; j < b.cols (); j++)
+              for (octave_idx_type k = 0; k < nc; k++)
+                for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                  retval(k,j) = b(ridx (i),j) / data (i);
+
+          if (calc_cond)
+            {
+              double dmax = 0., dmin = octave_Inf;
+              for (octave_idx_type i = 0; i < nm; i++)
+                {
+                  double tmp = fabs (data (i));
+                  if (tmp > dmax)
+                    dmax = tmp;
+                  if (tmp < dmin)
+                    dmin = tmp;
+                }
+              rcond = dmin / dmax;
+            }
+          else
+            rcond = 1.;
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseComplexMatrix
+SparseMatrix::dsolve (MatrixType &mattype, const SparseComplexMatrix& b,
+                     octave_idx_type& err, double& rcond,
+                     solve_singularity_handler, bool calc_cond) const
+{
+  SparseComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc < nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = SparseComplexMatrix (nc, b.cols ());
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Diagonal ||
+          typ == MatrixType::Permuted_Diagonal)
+        {
+          octave_idx_type b_nc = b.cols ();
+          octave_idx_type b_nz = b.nnz ();
+          retval = SparseComplexMatrix (nc, b_nc, b_nz);
+
+          retval.xcidx (0) = 0;
+          octave_idx_type ii = 0;
+          if (typ == MatrixType::Diagonal)
+            for (octave_idx_type j = 0; j < b.cols (); j++)
+              {
+                for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                  {
+                    if (b.ridx (i) >= nm)
+                      break;
+                    retval.xridx (ii) = b.ridx (i);
+                    retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
+                  }
+                retval.xcidx (j+1) = ii;
+              }
+          else
+            for (octave_idx_type j = 0; j < b.cols (); j++)
+              {
+                for (octave_idx_type l = 0; l < nc; l++)
+                  for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
+                    {
+                      bool found = false;
+                      octave_idx_type k;
+                      for (k = b.cidx (j); k < b.cidx (j+1); k++)
+                        if (ridx (i) == b.ridx (k))
+                          {
+                            found = true;
+                            break;
+                          }
+                      if (found)
+                        {
+                          retval.xridx (ii) = l;
+                          retval.xdata (ii++) = b.data (k) / data (i);
+                        }
+                    }
+                retval.xcidx (j+1) = ii;
+              }
+
+          if (calc_cond)
+            {
+              double dmax = 0., dmin = octave_Inf;
+              for (octave_idx_type i = 0; i < nm; i++)
+                {
+                  double tmp = fabs (data (i));
+                  if (tmp > dmax)
+                    dmax = tmp;
+                  if (tmp < dmin)
+                    dmin = tmp;
+                }
+              rcond = dmin / dmax;
+            }
+          else
+            rcond = 1.;
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+Matrix
+SparseMatrix::utsolve (MatrixType &mattype, const Matrix& b,
+                       octave_idx_type& err, double& rcond,
+                       solve_singularity_handler sing_handler,
+                       bool calc_cond) const
+{
+  Matrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc > nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = Matrix (nc, b.cols (), 0.0);
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Permuted_Upper ||
+          typ == MatrixType::Upper)
+        {
+          double anorm = 0.;
+          double ainvnorm = 0.;
+          octave_idx_type b_nc = b.cols ();
+          rcond = 1.;
+
+          if (calc_cond)
+            {
+              // Calculate the 1-norm of matrix for rcond calculation
+              for (octave_idx_type j = 0; j < nc; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          if (typ == MatrixType::Permuted_Upper)
+            {
+              retval.resize (nc, b_nc);
+              octave_idx_type *perm = mattype.triangular_perm ();
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    work[i] = b(i,j);
+                  for (octave_idx_type i = nr; i < nc; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                    {
+                      octave_idx_type kidx = perm[k];
+
+                      if (work[k] != 0.)
+                        {
+                          if (ridx (cidx (kidx+1)-1) != k ||
+                              data (cidx (kidx+1)-1) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          double tmp = work[k] / data (cidx (kidx+1)-1);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (kidx);
+                               i < cidx (kidx+1)-1; i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    retval.xelem (perm[i], j) = work[i];
+                }
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k >= 0; k--)
+                        {
+                          octave_idx_type iidx = perm[k];
+
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (iidx+1)-1);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (iidx);
+                                   i < cidx (iidx+1)-1; i++)
+                                {
+                                  octave_idx_type idx2 = ridx (i);
+                                  work[idx2] = work[idx2] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = 0; i < j+1; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+          else
+            {
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              retval.resize (nc, b_nc);
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    work[i] = b(i,j);
+                  for (octave_idx_type i = nr; i < nc; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                    {
+                      if (work[k] != 0.)
+                        {
+                          if (ridx (cidx (k+1)-1) != k ||
+                              data (cidx (k+1)-1) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          double tmp = work[k] / data (cidx (k+1)-1);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    retval.xelem (i, j) = work[i];
+                }
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k >= 0; k--)
+                        {
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (k+1)-1);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                                {
+                                  octave_idx_type iidx = ridx (i);
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = 0; i < j+1; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+
+        triangular_error:
+          if (err != 0)
+            {
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+
+          volatile double rcond_plus_one = rcond + 1.0;
+
+          if (rcond_plus_one == 1.0 || xisnan (rcond))
+            {
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseMatrix
+SparseMatrix::utsolve (MatrixType &mattype, const SparseMatrix& b,
+                       octave_idx_type& err, double& rcond,
+                       solve_singularity_handler sing_handler,
+                       bool calc_cond) const
+{
+  SparseMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc > nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = SparseMatrix (nc, b.cols ());
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Permuted_Upper ||
+          typ == MatrixType::Upper)
+        {
+          double anorm = 0.;
+          double ainvnorm = 0.;
+          rcond = 1.;
+
+          if (calc_cond)
+            {
+              // Calculate the 1-norm of matrix for rcond calculation
+              for (octave_idx_type j = 0; j < nc; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          octave_idx_type b_nc = b.cols ();
+          octave_idx_type b_nz = b.nnz ();
+          retval = SparseMatrix (nc, b_nc, b_nz);
+          retval.xcidx (0) = 0;
+          octave_idx_type ii = 0;
+          octave_idx_type x_nz = b_nz;
+
+          if (typ == MatrixType::Permuted_Upper)
+            {
+              octave_idx_type *perm = mattype.triangular_perm ();
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+
+              OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
+              for (octave_idx_type i = 0; i < nc; i++)
+                rperm[perm[i]] = i;
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                    work[b.ridx (i)] = b.data (i);
+
+                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                    {
+                      octave_idx_type kidx = perm[k];
+
+                      if (work[k] != 0.)
+                        {
+                          if (ridx (cidx (kidx+1)-1) != k ||
+                              data (cidx (kidx+1)-1) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          double tmp = work[k] / data (cidx (kidx+1)-1);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (kidx);
+                               i < cidx (kidx+1)-1; i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  // Count non-zeros in work vector and adjust space in
+                  // retval if needed
+                  octave_idx_type new_nnz = 0;
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (work[i] != 0.)
+                      new_nnz++;
+
+                  if (ii + new_nnz > x_nz)
+                    {
+                      // Resize the sparse matrix
+                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                      retval.change_capacity (sz);
+                      x_nz = sz;
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (work[rperm[i]] != 0.)
+                      {
+                        retval.xridx (ii) = i;
+                        retval.xdata (ii++) = work[rperm[i]];
+                      }
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k >= 0; k--)
+                        {
+                          octave_idx_type iidx = perm[k];
+
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (iidx+1)-1);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (iidx);
+                                   i < cidx (iidx+1)-1; i++)
+                                {
+                                  octave_idx_type idx2 = ridx (i);
+                                  work[idx2] = work[idx2] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = 0; i < j+1; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+          else
+            {
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                    work[b.ridx (i)] = b.data (i);
+
+                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                    {
+                      if (work[k] != 0.)
+                        {
+                          if (ridx (cidx (k+1)-1) != k ||
+                              data (cidx (k+1)-1) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          double tmp = work[k] / data (cidx (k+1)-1);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  // Count non-zeros in work vector and adjust space in
+                  // retval if needed
+                  octave_idx_type new_nnz = 0;
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (work[i] != 0.)
+                      new_nnz++;
+
+                  if (ii + new_nnz > x_nz)
+                    {
+                      // Resize the sparse matrix
+                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                      retval.change_capacity (sz);
+                      x_nz = sz;
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (work[i] != 0.)
+                      {
+                        retval.xridx (ii) = i;
+                        retval.xdata (ii++) = work[i];
+                      }
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k >= 0; k--)
+                        {
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (k+1)-1);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k);
+                                   i < cidx (k+1)-1; i++)
+                                {
+                                  octave_idx_type iidx = ridx (i);
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = 0; i < j+1; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+
+        triangular_error:
+          if (err != 0)
+            {
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+
+          volatile double rcond_plus_one = rcond + 1.0;
+
+          if (rcond_plus_one == 1.0 || xisnan (rcond))
+            {
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+  return retval;
+}
+
+ComplexMatrix
+SparseMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b,
+                       octave_idx_type& err, double& rcond,
+                       solve_singularity_handler sing_handler,
+                       bool calc_cond) const
+{
+  ComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc > nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Permuted_Upper ||
+          typ == MatrixType::Upper)
+        {
+          double anorm = 0.;
+          double ainvnorm = 0.;
+          octave_idx_type b_nc = b.cols ();
+          rcond = 1.;
+
+          if (calc_cond)
+            {
+              // Calculate the 1-norm of matrix for rcond calculation
+              for (octave_idx_type j = 0; j < nc; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          if (typ == MatrixType::Permuted_Upper)
+            {
+              retval.resize (nc, b_nc);
+              octave_idx_type *perm = mattype.triangular_perm ();
+              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    cwork[i] = b(i,j);
+                  for (octave_idx_type i = nr; i < nc; i++)
+                    cwork[i] = 0.;
+
+                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                    {
+                      octave_idx_type kidx = perm[k];
+
+                      if (cwork[k] != 0.)
+                        {
+                          if (ridx (cidx (kidx+1)-1) != k ||
+                              data (cidx (kidx+1)-1) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
+                          cwork[k] = tmp;
+                          for (octave_idx_type i = cidx (kidx);
+                               i < cidx (kidx+1)-1; i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    retval.xelem (perm[i], j) = cwork[i];
+                }
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  OCTAVE_LOCAL_BUFFER (double, work, nm);
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k >= 0; k--)
+                        {
+                          octave_idx_type iidx = perm[k];
+
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (iidx+1)-1);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (iidx);
+                                   i < cidx (iidx+1)-1; i++)
+                                {
+                                  octave_idx_type idx2 = ridx (i);
+                                  work[idx2] = work[idx2] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = 0; i < j+1; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+          else
+            {
+              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+              retval.resize (nc, b_nc);
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    cwork[i] = b(i,j);
+                  for (octave_idx_type i = nr; i < nc; i++)
+                    cwork[i] = 0.;
+
+                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                    {
+                      if (cwork[k] != 0.)
+                        {
+                          if (ridx (cidx (k+1)-1) != k ||
+                              data (cidx (k+1)-1) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          Complex tmp = cwork[k] / data (cidx (k+1)-1);
+                          cwork[k] = tmp;
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              cwork[iidx] = cwork[iidx] - tmp  * data (i);
+                            }
+                        }
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    retval.xelem (i, j) = cwork[i];
+                }
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  OCTAVE_LOCAL_BUFFER (double, work, nm);
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k >= 0; k--)
+                        {
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (k+1)-1);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k);
+                                   i < cidx (k+1)-1; i++)
+                                {
+                                  octave_idx_type iidx = ridx (i);
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = 0; i < j+1; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+
+        triangular_error:
+          if (err != 0)
+            {
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+
+          volatile double rcond_plus_one = rcond + 1.0;
+
+          if (rcond_plus_one == 1.0 || xisnan (rcond))
+            {
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseComplexMatrix
+SparseMatrix::utsolve (MatrixType &mattype, const SparseComplexMatrix& b,
+                       octave_idx_type& err, double& rcond,
+                       solve_singularity_handler sing_handler,
+                       bool calc_cond) const
+{
+  SparseComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc > nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = SparseComplexMatrix (nc, b.cols ());
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Permuted_Upper ||
+          typ == MatrixType::Upper)
+        {
+          double anorm = 0.;
+          double ainvnorm = 0.;
+          rcond = 1.;
+
+          if (calc_cond)
+            {
+              // Calculate the 1-norm of matrix for rcond calculation
+              for (octave_idx_type j = 0; j < nc; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          octave_idx_type b_nc = b.cols ();
+          octave_idx_type b_nz = b.nnz ();
+          retval = SparseComplexMatrix (nc, b_nc, b_nz);
+          retval.xcidx (0) = 0;
+          octave_idx_type ii = 0;
+          octave_idx_type x_nz = b_nz;
+
+          if (typ == MatrixType::Permuted_Upper)
+            {
+              octave_idx_type *perm = mattype.triangular_perm ();
+              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+
+              OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
+              for (octave_idx_type i = 0; i < nc; i++)
+                rperm[perm[i]] = i;
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    cwork[i] = 0.;
+                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                    cwork[b.ridx (i)] = b.data (i);
+
+                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                    {
+                      octave_idx_type kidx = perm[k];
+
+                      if (cwork[k] != 0.)
+                        {
+                          if (ridx (cidx (kidx+1)-1) != k ||
+                              data (cidx (kidx+1)-1) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
+                          cwork[k] = tmp;
+                          for (octave_idx_type i = cidx (kidx);
+                               i < cidx (kidx+1)-1; i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  // Count non-zeros in work vector and adjust space in
+                  // retval if needed
+                  octave_idx_type new_nnz = 0;
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (cwork[i] != 0.)
+                      new_nnz++;
+
+                  if (ii + new_nnz > x_nz)
+                    {
+                      // Resize the sparse matrix
+                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                      retval.change_capacity (sz);
+                      x_nz = sz;
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (cwork[rperm[i]] != 0.)
+                      {
+                        retval.xridx (ii) = i;
+                        retval.xdata (ii++) = cwork[rperm[i]];
+                      }
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  OCTAVE_LOCAL_BUFFER (double, work, nm);
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k >= 0; k--)
+                        {
+                          octave_idx_type iidx = perm[k];
+
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (iidx+1)-1);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (iidx);
+                                   i < cidx (iidx+1)-1; i++)
+                                {
+                                  octave_idx_type idx2 = ridx (i);
+                                  work[idx2] = work[idx2] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = 0; i < j+1; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+          else
+            {
+              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    cwork[i] = 0.;
+                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                    cwork[b.ridx (i)] = b.data (i);
+
+                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                    {
+                      if (cwork[k] != 0.)
+                        {
+                          if (ridx (cidx (k+1)-1) != k ||
+                              data (cidx (k+1)-1) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          Complex tmp = cwork[k] / data (cidx (k+1)-1);
+                          cwork[k] = tmp;
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  // Count non-zeros in work vector and adjust space in
+                  // retval if needed
+                  octave_idx_type new_nnz = 0;
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (cwork[i] != 0.)
+                      new_nnz++;
+
+                  if (ii + new_nnz > x_nz)
+                    {
+                      // Resize the sparse matrix
+                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                      retval.change_capacity (sz);
+                      x_nz = sz;
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (cwork[i] != 0.)
+                      {
+                        retval.xridx (ii) = i;
+                        retval.xdata (ii++) = cwork[i];
+                      }
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  OCTAVE_LOCAL_BUFFER (double, work, nm);
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k >= 0; k--)
+                        {
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (k+1)-1);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k);
+                                   i < cidx (k+1)-1; i++)
+                                {
+                                  octave_idx_type iidx = ridx (i);
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = 0; i < j+1; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+
+        triangular_error:
+          if (err != 0)
+            {
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+
+          volatile double rcond_plus_one = rcond + 1.0;
+
+          if (rcond_plus_one == 1.0 || xisnan (rcond))
+            {
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+Matrix
+SparseMatrix::ltsolve (MatrixType &mattype, const Matrix& b,
+                       octave_idx_type& err, double& rcond,
+                       solve_singularity_handler sing_handler,
+                       bool calc_cond) const
+{
+  Matrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc > nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = Matrix (nc, b.cols (), 0.0);
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Permuted_Lower ||
+          typ == MatrixType::Lower)
+        {
+          double anorm = 0.;
+          double ainvnorm = 0.;
+          octave_idx_type b_nc = b.cols ();
+          rcond = 1.;
+
+          if (calc_cond)
+            {
+              // Calculate the 1-norm of matrix for rcond calculation
+              for (octave_idx_type j = 0; j < nc; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          if (typ == MatrixType::Permuted_Lower)
+            {
+              retval.resize (nc, b_nc);
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              octave_idx_type *perm = mattype.triangular_perm ();
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  if (nc > nr)
+                    for (octave_idx_type i = 0; i < nm; i++)
+                      work[i] = 0.;
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    work[perm[i]] = b(i,j);
+
+                  for (octave_idx_type k = 0; k < nc; k++)
+                    {
+                      if (work[k] != 0.)
+                        {
+                          octave_idx_type minr = nr;
+                          octave_idx_type mini = 0;
+
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                            if (perm[ridx (i)] < minr)
+                              {
+                                minr = perm[ridx (i)];
+                                mini = i;
+                              }
+
+                          if (minr != k || data (mini) == 0)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          double tmp = work[k] / data (mini);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                            {
+                              if (i == mini)
+                                continue;
+
+                              octave_idx_type iidx = perm[ridx (i)];
+                              work[iidx] = work[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    retval(i, j) = work[i];
+                }
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = 0; k < nc; k++)
+                        {
+                          if (work[k] != 0.)
+                            {
+                              octave_idx_type minr = nr;
+                              octave_idx_type mini = 0;
+
+                              for (octave_idx_type i = cidx (k);
+                                   i < cidx (k+1); i++)
+                                if (perm[ridx (i)] < minr)
+                                  {
+                                    minr = perm[ridx (i)];
+                                    mini = i;
+                                  }
+
+                              double tmp = work[k] / data (mini);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k);
+                                   i < cidx (k+1); i++)
+                                {
+                                  if (i == mini)
+                                    continue;
+
+                                  octave_idx_type iidx = perm[ridx (i)];
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+
+                      double atmp = 0;
+                      for (octave_idx_type i = j; i < nc; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+          else
+            {
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              retval.resize (nc, b_nc, 0.);
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    work[i] = b(i,j);
+                  for (octave_idx_type i = nr; i < nc; i++)
+                    work[i] = 0.;
+                  for (octave_idx_type k = 0; k < nc; k++)
+                    {
+                      if (work[k] != 0.)
+                        {
+                          if (ridx (cidx (k)) != k ||
+                              data (cidx (k)) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          double tmp = work[k] / data (cidx (k));
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k)+1;
+                               i < cidx (k+1); i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    retval.xelem (i, j) = work[i];
+                }
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k < nc; k++)
+                        {
+
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (k));
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k)+1;
+                                   i < cidx (k+1); i++)
+                                {
+                                  octave_idx_type iidx = ridx (i);
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = j; i < nc; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+
+        triangular_error:
+          if (err != 0)
+            {
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+
+          volatile double rcond_plus_one = rcond + 1.0;
+
+          if (rcond_plus_one == 1.0 || xisnan (rcond))
+            {
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseMatrix
+SparseMatrix::ltsolve (MatrixType &mattype, const SparseMatrix& b,
+                       octave_idx_type& err, double& rcond,
+                       solve_singularity_handler sing_handler,
+                       bool calc_cond) const
+{
+  SparseMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc > nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = SparseMatrix (nc, b.cols ());
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Permuted_Lower ||
+          typ == MatrixType::Lower)
+        {
+          double anorm = 0.;
+          double ainvnorm = 0.;
+          rcond = 1.;
+
+          if (calc_cond)
+            {
+              // Calculate the 1-norm of matrix for rcond calculation
+              for (octave_idx_type j = 0; j < nc; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          octave_idx_type b_nc = b.cols ();
+          octave_idx_type b_nz = b.nnz ();
+          retval = SparseMatrix (nc, b_nc, b_nz);
+          retval.xcidx (0) = 0;
+          octave_idx_type ii = 0;
+          octave_idx_type x_nz = b_nz;
+
+          if (typ == MatrixType::Permuted_Lower)
+            {
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              octave_idx_type *perm = mattype.triangular_perm ();
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                    work[perm[b.ridx (i)]] = b.data (i);
+
+                  for (octave_idx_type k = 0; k < nc; k++)
+                    {
+                      if (work[k] != 0.)
+                        {
+                          octave_idx_type minr = nr;
+                          octave_idx_type mini = 0;
+
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                            if (perm[ridx (i)] < minr)
+                              {
+                                minr = perm[ridx (i)];
+                                mini = i;
+                              }
+
+                          if (minr != k || data (mini) == 0)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          double tmp = work[k] / data (mini);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                            {
+                              if (i == mini)
+                                continue;
+
+                              octave_idx_type iidx = perm[ridx (i)];
+                              work[iidx] = work[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  // Count non-zeros in work vector and adjust space in
+                  // retval if needed
+                  octave_idx_type new_nnz = 0;
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (work[i] != 0.)
+                      new_nnz++;
+
+                  if (ii + new_nnz > x_nz)
+                    {
+                      // Resize the sparse matrix
+                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                      retval.change_capacity (sz);
+                      x_nz = sz;
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (work[i] != 0.)
+                      {
+                        retval.xridx (ii) = i;
+                        retval.xdata (ii++) = work[i];
+                      }
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = 0; k < nc; k++)
+                        {
+                          if (work[k] != 0.)
+                            {
+                              octave_idx_type minr = nr;
+                              octave_idx_type mini = 0;
+
+                              for (octave_idx_type i = cidx (k);
+                                   i < cidx (k+1); i++)
+                                if (perm[ridx (i)] < minr)
+                                  {
+                                    minr = perm[ridx (i)];
+                                    mini = i;
+                                  }
+
+                              double tmp = work[k] / data (mini);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k);
+                                   i < cidx (k+1); i++)
+                                {
+                                  if (i == mini)
+                                    continue;
+
+                                  octave_idx_type iidx = perm[ridx (i)];
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+
+                      double atmp = 0;
+                      for (octave_idx_type i = j; i < nr; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+          else
+            {
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                    work[b.ridx (i)] = b.data (i);
+
+                  for (octave_idx_type k = 0; k < nc; k++)
+                    {
+                      if (work[k] != 0.)
+                        {
+                          if (ridx (cidx (k)) != k ||
+                              data (cidx (k)) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          double tmp = work[k] / data (cidx (k));
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  // Count non-zeros in work vector and adjust space in
+                  // retval if needed
+                  octave_idx_type new_nnz = 0;
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (work[i] != 0.)
+                      new_nnz++;
+
+                  if (ii + new_nnz > x_nz)
+                    {
+                      // Resize the sparse matrix
+                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                      retval.change_capacity (sz);
+                      x_nz = sz;
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (work[i] != 0.)
+                      {
+                        retval.xridx (ii) = i;
+                        retval.xdata (ii++) = work[i];
+                      }
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k < nc; k++)
+                        {
+
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (k));
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k)+1;
+                                   i < cidx (k+1); i++)
+                                {
+                                  octave_idx_type iidx = ridx (i);
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = j; i < nc; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+
+        triangular_error:
+          if (err != 0)
+            {
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+
+          volatile double rcond_plus_one = rcond + 1.0;
+
+          if (rcond_plus_one == 1.0 || xisnan (rcond))
+            {
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+ComplexMatrix
+SparseMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b,
+                       octave_idx_type& err, double& rcond,
+                       solve_singularity_handler sing_handler,
+                       bool calc_cond) const
+{
+  ComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc > nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Permuted_Lower ||
+          typ == MatrixType::Lower)
+        {
+          double anorm = 0.;
+          double ainvnorm = 0.;
+          octave_idx_type b_nc = b.cols ();
+          rcond = 1.;
+
+          if (calc_cond)
+            {
+              // Calculate the 1-norm of matrix for rcond calculation
+              for (octave_idx_type j = 0; j < nc; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          if (typ == MatrixType::Permuted_Lower)
+            {
+              retval.resize (nc, b_nc);
+              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+              octave_idx_type *perm = mattype.triangular_perm ();
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    cwork[i] = 0.;
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    cwork[perm[i]] = b(i,j);
+
+                  for (octave_idx_type k = 0; k < nc; k++)
+                    {
+                      if (cwork[k] != 0.)
+                        {
+                          octave_idx_type minr = nr;
+                          octave_idx_type mini = 0;
+
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                            if (perm[ridx (i)] < minr)
+                              {
+                                minr = perm[ridx (i)];
+                                mini = i;
+                              }
+
+                          if (minr != k || data (mini) == 0)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          Complex tmp = cwork[k] / data (mini);
+                          cwork[k] = tmp;
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                            {
+                              if (i == mini)
+                                continue;
+
+                              octave_idx_type iidx = perm[ridx (i)];
+                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    retval(i, j) = cwork[i];
+                }
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  OCTAVE_LOCAL_BUFFER (double, work, nm);
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = 0; k < nc; k++)
+                        {
+                          if (work[k] != 0.)
+                            {
+                              octave_idx_type minr = nr;
+                              octave_idx_type mini = 0;
+
+                              for (octave_idx_type i = cidx (k);
+                                   i < cidx (k+1); i++)
+                                if (perm[ridx (i)] < minr)
+                                  {
+                                    minr = perm[ridx (i)];
+                                    mini = i;
+                                  }
+
+                              double tmp = work[k] / data (mini);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k);
+                                   i < cidx (k+1); i++)
+                                {
+                                  if (i == mini)
+                                    continue;
+
+                                  octave_idx_type iidx = perm[ridx (i)];
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+
+                      double atmp = 0;
+                      for (octave_idx_type i = j; i < nc; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+          else
+            {
+              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+              retval.resize (nc, b_nc, 0.);
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    cwork[i] = b(i,j);
+                  for (octave_idx_type i = nr; i < nc; i++)
+                    cwork[i] = 0.;
+
+                  for (octave_idx_type k = 0; k < nc; k++)
+                    {
+                      if (cwork[k] != 0.)
+                        {
+                          if (ridx (cidx (k)) != k ||
+                              data (cidx (k)) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          Complex tmp = cwork[k] / data (cidx (k));
+                          cwork[k] = tmp;
+                          for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    retval.xelem (i, j) = cwork[i];
+                }
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  OCTAVE_LOCAL_BUFFER (double, work, nm);
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k < nc; k++)
+                        {
+
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (k));
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k)+1;
+                                   i < cidx (k+1); i++)
+                                {
+                                  octave_idx_type iidx = ridx (i);
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = j; i < nc; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+
+        triangular_error:
+          if (err != 0)
+            {
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+
+          volatile double rcond_plus_one = rcond + 1.0;
+
+          if (rcond_plus_one == 1.0 || xisnan (rcond))
+            {
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseComplexMatrix
+SparseMatrix::ltsolve (MatrixType &mattype, const SparseComplexMatrix& b,
+                       octave_idx_type& err, double& rcond,
+                       solve_singularity_handler sing_handler,
+                       bool calc_cond) const
+{
+  SparseComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nm = (nc > nr ? nc : nr);
+  err = 0;
+
+  if (nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || nc == 0 || b.cols () == 0)
+    retval = SparseComplexMatrix (nc, b.cols ());
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Permuted_Lower ||
+          typ == MatrixType::Lower)
+        {
+          double anorm = 0.;
+          double ainvnorm = 0.;
+          rcond = 1.;
+
+          if (calc_cond)
+            {
+              // Calculate the 1-norm of matrix for rcond calculation
+              for (octave_idx_type j = 0; j < nc; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          octave_idx_type b_nc = b.cols ();
+          octave_idx_type b_nz = b.nnz ();
+          retval = SparseComplexMatrix (nc, b_nc, b_nz);
+          retval.xcidx (0) = 0;
+          octave_idx_type ii = 0;
+          octave_idx_type x_nz = b_nz;
+
+          if (typ == MatrixType::Permuted_Lower)
+            {
+              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+              octave_idx_type *perm = mattype.triangular_perm ();
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    cwork[i] = 0.;
+                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                    cwork[perm[b.ridx (i)]] = b.data (i);
+
+                  for (octave_idx_type k = 0; k < nc; k++)
+                    {
+                      if (cwork[k] != 0.)
+                        {
+                          octave_idx_type minr = nr;
+                          octave_idx_type mini = 0;
+
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                            if (perm[ridx (i)] < minr)
+                              {
+                                minr = perm[ridx (i)];
+                                mini = i;
+                              }
+
+                          if (minr != k || data (mini) == 0)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          Complex tmp = cwork[k] / data (mini);
+                          cwork[k] = tmp;
+                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                            {
+                              if (i == mini)
+                                continue;
+
+                              octave_idx_type iidx = perm[ridx (i)];
+                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  // Count non-zeros in work vector and adjust space in
+                  // retval if needed
+                  octave_idx_type new_nnz = 0;
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (cwork[i] != 0.)
+                      new_nnz++;
+
+                  if (ii + new_nnz > x_nz)
+                    {
+                      // Resize the sparse matrix
+                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                      retval.change_capacity (sz);
+                      x_nz = sz;
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (cwork[i] != 0.)
+                      {
+                        retval.xridx (ii) = i;
+                        retval.xdata (ii++) = cwork[i];
+                      }
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  OCTAVE_LOCAL_BUFFER (double, work, nm);
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = 0; k < nc; k++)
+                        {
+                          if (work[k] != 0.)
+                            {
+                              octave_idx_type minr = nr;
+                              octave_idx_type mini = 0;
+
+                              for (octave_idx_type i = cidx (k);
+                                   i < cidx (k+1); i++)
+                                if (perm[ridx (i)] < minr)
+                                  {
+                                    minr = perm[ridx (i)];
+                                    mini = i;
+                                  }
+
+                              double tmp = work[k] / data (mini);
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k);
+                                   i < cidx (k+1); i++)
+                                {
+                                  if (i == mini)
+                                    continue;
+
+                                  octave_idx_type iidx = perm[ridx (i)];
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+
+                      double atmp = 0;
+                      for (octave_idx_type i = j; i < nc; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+          else
+            {
+              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    cwork[i] = 0.;
+                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                    cwork[b.ridx (i)] = b.data (i);
+
+                  for (octave_idx_type k = 0; k < nc; k++)
+                    {
+                      if (cwork[k] != 0.)
+                        {
+                          if (ridx (cidx (k)) != k ||
+                              data (cidx (k)) == 0.)
+                            {
+                              err = -2;
+                              goto triangular_error;
+                            }
+
+                          Complex tmp = cwork[k] / data (cidx (k));
+                          cwork[k] = tmp;
+                          for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                            {
+                              octave_idx_type iidx = ridx (i);
+                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                            }
+                        }
+                    }
+
+                  // Count non-zeros in work vector and adjust space in
+                  // retval if needed
+                  octave_idx_type new_nnz = 0;
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (cwork[i] != 0.)
+                      new_nnz++;
+
+                  if (ii + new_nnz > x_nz)
+                    {
+                      // Resize the sparse matrix
+                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                      retval.change_capacity (sz);
+                      x_nz = sz;
+                    }
+
+                  for (octave_idx_type i = 0; i < nc; i++)
+                    if (cwork[i] != 0.)
+                      {
+                        retval.xridx (ii) = i;
+                        retval.xdata (ii++) = cwork[i];
+                      }
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+
+              if (calc_cond)
+                {
+                  // Calculation of 1-norm of inv(*this)
+                  OCTAVE_LOCAL_BUFFER (double, work, nm);
+                  for (octave_idx_type i = 0; i < nm; i++)
+                    work[i] = 0.;
+
+                  for (octave_idx_type j = 0; j < nr; j++)
+                    {
+                      work[j] = 1.;
+
+                      for (octave_idx_type k = j; k < nc; k++)
+                        {
+
+                          if (work[k] != 0.)
+                            {
+                              double tmp = work[k] / data (cidx (k));
+                              work[k] = tmp;
+                              for (octave_idx_type i = cidx (k)+1;
+                                   i < cidx (k+1); i++)
+                                {
+                                  octave_idx_type iidx = ridx (i);
+                                  work[iidx] = work[iidx] - tmp * data (i);
+                                }
+                            }
+                        }
+                      double atmp = 0;
+                      for (octave_idx_type i = j; i < nc; i++)
+                        {
+                          atmp += fabs (work[i]);
+                          work[i] = 0.;
+                        }
+                      if (atmp > ainvnorm)
+                        ainvnorm = atmp;
+                    }
+                  rcond = 1. / ainvnorm / anorm;
+                }
+            }
+
+        triangular_error:
+          if (err != 0)
+            {
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+
+          volatile double rcond_plus_one = rcond + 1.0;
+
+          if (rcond_plus_one == 1.0 || xisnan (rcond))
+            {
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision, rcond = %g",
+                   rcond);
+            }
+        }
+      else
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+Matrix
+SparseMatrix::trisolve (MatrixType &mattype, const Matrix& b,
+                        octave_idx_type& err, double& rcond,
+                        solve_singularity_handler sing_handler,
+                        bool calc_cond) const
+{
+  Matrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = Matrix (nc, b.cols (), 0.0);
+  else if (calc_cond)
+    (*current_liboctave_error_handler)
+      ("calculation of condition number not implemented");
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      volatile int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Tridiagonal_Hermitian)
+        {
+          OCTAVE_LOCAL_BUFFER (double, D, nr);
+          OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
+
+          if (mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < nc-1; j++)
+                {
+                  D[j] = data (ii++);
+                  DL[j] = data (ii);
+                  ii += 2;
+                }
+              D[nc-1] = data (ii);
+            }
+          else
+            {
+              D[0] = 0.;
+              for (octave_idx_type i = 0; i < nr - 1; i++)
+                {
+                  D[i+1] = 0.;
+                  DL[i] = 0.;
+                }
+
+              for (octave_idx_type j = 0; j < nc; j++)
+                for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                  {
+                    if (ridx (i) == j)
+                      D[j] = data (i);
+                    else if (ridx (i) == j + 1)
+                      DL[j] = data (i);
+                  }
+            }
+
+          octave_idx_type b_nc = b.cols ();
+          retval = b;
+          double *result = retval.fortran_vec ();
+
+          F77_XFCN (dptsv, DPTSV, (nr, b_nc, D, DL, result,
+                                   b.rows (), err));
+
+          if (err != 0)
+            {
+              err = 0;
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Tridiagonal;
+            }
+          else
+            rcond = 1.;
+        }
+
+      if (typ == MatrixType::Tridiagonal)
+        {
+          OCTAVE_LOCAL_BUFFER (double, DU, nr - 1);
+          OCTAVE_LOCAL_BUFFER (double, D, nr);
+          OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
+
+          if (mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < nc-1; j++)
+                {
+                  D[j] = data (ii++);
+                  DL[j] = data (ii++);
+                  DU[j] = data (ii++);
+                }
+              D[nc-1] = data (ii);
+            }
+          else
+            {
+              D[0] = 0.;
+              for (octave_idx_type i = 0; i < nr - 1; i++)
+                {
+                  D[i+1] = 0.;
+                  DL[i] = 0.;
+                  DU[i] = 0.;
+                }
+
+              for (octave_idx_type j = 0; j < nc; j++)
+                for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                  {
+                    if (ridx (i) == j)
+                      D[j] = data (i);
+                    else if (ridx (i) == j + 1)
+                      DL[j] = data (i);
+                    else if (ridx (i) == j - 1)
+                      DU[j-1] = data (i);
+                  }
+            }
+
+          octave_idx_type b_nc = b.cols ();
+          retval = b;
+          double *result = retval.fortran_vec ();
+
+          F77_XFCN (dgtsv, DGTSV, (nr, b_nc, DL, D, DU, result,
+                                   b.rows (), err));
+
+          if (err != 0)
+            {
+              rcond = 0.;
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision");
+
+            }
+          else
+            rcond = 1.;
+        }
+      else if (typ != MatrixType::Tridiagonal_Hermitian)
+               (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseMatrix
+SparseMatrix::trisolve (MatrixType &mattype, const SparseMatrix& b,
+                        octave_idx_type& err, double& rcond,
+                        solve_singularity_handler sing_handler,
+                        bool calc_cond) const
+{
+  SparseMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = SparseMatrix (nc, b.cols ());
+  else if (calc_cond)
+    (*current_liboctave_error_handler)
+      ("calculation of condition number not implemented");
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      // Note can't treat symmetric case as there is no dpttrf function
+      if (typ == MatrixType::Tridiagonal ||
+          typ == MatrixType::Tridiagonal_Hermitian)
+        {
+          OCTAVE_LOCAL_BUFFER (double, DU2, nr - 2);
+          OCTAVE_LOCAL_BUFFER (double, DU, nr - 1);
+          OCTAVE_LOCAL_BUFFER (double, D, nr);
+          OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
+          Array<octave_idx_type> ipvt (dim_vector (nr, 1));
+          octave_idx_type *pipvt = ipvt.fortran_vec ();
+
+          if (mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < nc-1; j++)
+                {
+                  D[j] = data (ii++);
+                  DL[j] = data (ii++);
+                  DU[j] = data (ii++);
+                }
+              D[nc-1] = data (ii);
+            }
+          else
+            {
+              D[0] = 0.;
+              for (octave_idx_type i = 0; i < nr - 1; i++)
+                {
+                  D[i+1] = 0.;
+                  DL[i] = 0.;
+                  DU[i] = 0.;
+                }
+
+              for (octave_idx_type j = 0; j < nc; j++)
+                for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                  {
+                    if (ridx (i) == j)
+                      D[j] = data (i);
+                    else if (ridx (i) == j + 1)
+                      DL[j] = data (i);
+                    else if (ridx (i) == j - 1)
+                      DU[j-1] = data (i);
+                  }
+            }
+
+          F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
+
+          if (err != 0)
+            {
+              rcond = 0.0;
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision");
+
+            }
+          else
+            {
+              rcond = 1.0;
+              char job = 'N';
+              volatile octave_idx_type x_nz = b.nnz ();
+              octave_idx_type b_nc = b.cols ();
+              retval = SparseMatrix (nr, b_nc, x_nz);
+              retval.xcidx (0) = 0;
+              volatile octave_idx_type ii = 0;
+
+              OCTAVE_LOCAL_BUFFER (double, work, nr);
+
+              for (volatile octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    work[i] = 0.;
+                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                    work[b.ridx (i)] = b.data (i);
+
+                  F77_XFCN (dgttrs, DGTTRS,
+                            (F77_CONST_CHAR_ARG2 (&job, 1),
+                             nr, 1, DL, D, DU, DU2, pipvt,
+                             work, b.rows (), err
+                             F77_CHAR_ARG_LEN (1)));
+
+                  // Count non-zeros in work vector and adjust
+                  // space in retval if needed
+                  octave_idx_type new_nnz = 0;
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    if (work[i] != 0.)
+                      new_nnz++;
+
+                  if (ii + new_nnz > x_nz)
+                    {
+                      // Resize the sparse matrix
+                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                      retval.change_capacity (sz);
+                      x_nz = sz;
+                    }
+
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    if (work[i] != 0.)
+                      {
+                        retval.xridx (ii) = i;
+                        retval.xdata (ii++) = work[i];
+                      }
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+            }
+        }
+      else if (typ != MatrixType::Tridiagonal_Hermitian)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+ComplexMatrix
+SparseMatrix::trisolve (MatrixType &mattype, const ComplexMatrix& b,
+                        octave_idx_type& err, double& rcond,
+                        solve_singularity_handler sing_handler,
+                        bool calc_cond) const
+{
+  ComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
+  else if (calc_cond)
+    (*current_liboctave_error_handler)
+      ("calculation of condition number not implemented");
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      volatile int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Tridiagonal_Hermitian)
+        {
+          OCTAVE_LOCAL_BUFFER (double, D, nr);
+          OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
+
+          if (mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < nc-1; j++)
+                {
+                  D[j] = data (ii++);
+                  DL[j] = data (ii);
+                  ii += 2;
+                }
+              D[nc-1] = data (ii);
+            }
+          else
+            {
+              D[0] = 0.;
+              for (octave_idx_type i = 0; i < nr - 1; i++)
+                {
+                  D[i+1] = 0.;
+                  DL[i] = 0.;
+                }
+
+              for (octave_idx_type j = 0; j < nc; j++)
+                for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                  {
+                    if (ridx (i) == j)
+                      D[j] = data (i);
+                    else if (ridx (i) == j + 1)
+                      DL[j] = data (i);
+                  }
+            }
+
+          octave_idx_type b_nr = b.rows ();
+          octave_idx_type b_nc = b.cols ();
+          rcond = 1.;
+
+          retval = b;
+          Complex *result = retval.fortran_vec ();
+
+          F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result,
+                                   b_nr, err));
+
+          if (err != 0)
+            {
+              err = 0;
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Tridiagonal;
+            }
+        }
+
+      if (typ == MatrixType::Tridiagonal)
+        {
+          OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
+          OCTAVE_LOCAL_BUFFER (Complex, D, nr);
+          OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
+
+          if (mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < nc-1; j++)
+                {
+                  D[j] = data (ii++);
+                  DL[j] = data (ii++);
+                  DU[j] = data (ii++);
+                }
+              D[nc-1] = data (ii);
+            }
+          else
+            {
+              D[0] = 0.;
+              for (octave_idx_type i = 0; i < nr - 1; i++)
+                {
+                  D[i+1] = 0.;
+                  DL[i] = 0.;
+                  DU[i] = 0.;
+                }
+
+              for (octave_idx_type j = 0; j < nc; j++)
+                for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                  {
+                    if (ridx (i) == j)
+                      D[j] = data (i);
+                    else if (ridx (i) == j + 1)
+                      DL[j] = data (i);
+                    else if (ridx (i) == j - 1)
+                      DU[j-1] = data (i);
+                  }
+            }
+
+          octave_idx_type b_nr = b.rows ();
+          octave_idx_type b_nc = b.cols ();
+          rcond = 1.;
+
+          retval = b;
+          Complex *result = retval.fortran_vec ();
+
+          F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result,
+                                   b_nr, err));
+
+          if (err != 0)
+            {
+              rcond = 0.;
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision");
+            }
+        }
+      else if (typ != MatrixType::Tridiagonal_Hermitian)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseComplexMatrix
+SparseMatrix::trisolve (MatrixType &mattype, const SparseComplexMatrix& b,
+                        octave_idx_type& err, double& rcond,
+                        solve_singularity_handler sing_handler,
+                        bool calc_cond) const
+{
+  SparseComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = SparseComplexMatrix (nc, b.cols ());
+  else if (calc_cond)
+    (*current_liboctave_error_handler)
+      ("calculation of condition number not implemented");
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      int typ = mattype.type ();
+      mattype.info ();
+
+      // Note can't treat symmetric case as there is no dpttrf function
+      if (typ == MatrixType::Tridiagonal ||
+          typ == MatrixType::Tridiagonal_Hermitian)
+        {
+          OCTAVE_LOCAL_BUFFER (double, DU2, nr - 2);
+          OCTAVE_LOCAL_BUFFER (double, DU, nr - 1);
+          OCTAVE_LOCAL_BUFFER (double, D, nr);
+          OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
+          Array<octave_idx_type> ipvt (dim_vector (nr, 1));
+          octave_idx_type *pipvt = ipvt.fortran_vec ();
+
+          if (mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < nc-1; j++)
+                {
+                  D[j] = data (ii++);
+                  DL[j] = data (ii++);
+                  DU[j] = data (ii++);
+                }
+              D[nc-1] = data (ii);
+            }
+          else
+            {
+              D[0] = 0.;
+              for (octave_idx_type i = 0; i < nr - 1; i++)
+                {
+                  D[i+1] = 0.;
+                  DL[i] = 0.;
+                  DU[i] = 0.;
+                }
+
+              for (octave_idx_type j = 0; j < nc; j++)
+                for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                  {
+                    if (ridx (i) == j)
+                      D[j] = data (i);
+                    else if (ridx (i) == j + 1)
+                      DL[j] = data (i);
+                    else if (ridx (i) == j - 1)
+                      DU[j-1] = data (i);
+                  }
+            }
+
+          F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
+
+          if (err != 0)
+            {
+              rcond = 0.0;
+              err = -2;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision");
+            }
+          else
+            {
+              rcond = 1.;
+              char job = 'N';
+              octave_idx_type b_nr = b.rows ();
+              octave_idx_type b_nc = b.cols ();
+              OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
+              OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
+
+              // Take a first guess that the number of non-zero terms
+              // will be as many as in b
+              volatile octave_idx_type x_nz = b.nnz ();
+              volatile octave_idx_type ii = 0;
+              retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
+
+              retval.xcidx (0) = 0;
+              for (volatile octave_idx_type j = 0; j < b_nc; j++)
+                {
+
+                  for (octave_idx_type i = 0; i < b_nr; i++)
+                    {
+                      Complex c = b (i,j);
+                      Bx[i] = std::real (c);
+                      Bz[i] = std::imag (c);
+                    }
+
+                  F77_XFCN (dgttrs, DGTTRS,
+                            (F77_CONST_CHAR_ARG2 (&job, 1),
+                             nr, 1, DL, D, DU, DU2, pipvt,
+                             Bx, b_nr, err
+                             F77_CHAR_ARG_LEN (1)));
+
+                  if (err != 0)
+                    {
+                      (*current_liboctave_error_handler)
+                        ("SparseMatrix::solve solve failed");
+
+                      err = -1;
+                      break;
+                    }
+
+                  F77_XFCN (dgttrs, DGTTRS,
+                            (F77_CONST_CHAR_ARG2 (&job, 1),
+                             nr, 1, DL, D, DU, DU2, pipvt,
+                             Bz, b_nr, err
+                             F77_CHAR_ARG_LEN (1)));
+
+                  if (err != 0)
+                    {
+                      (*current_liboctave_error_handler)
+                        ("SparseMatrix::solve solve failed");
+
+                      err = -1;
+                      break;
+                    }
+
+                  // Count non-zeros in work vector and adjust
+                  // space in retval if needed
+                  octave_idx_type new_nnz = 0;
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    if (Bx[i] != 0. || Bz[i] != 0.)
+                      new_nnz++;
+
+                  if (ii + new_nnz > x_nz)
+                    {
+                      // Resize the sparse matrix
+                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                      retval.change_capacity (sz);
+                      x_nz = sz;
+                    }
+
+                  for (octave_idx_type i = 0; i < nr; i++)
+                    if (Bx[i] != 0. || Bz[i] != 0.)
+                      {
+                        retval.xridx (ii) = i;
+                        retval.xdata (ii++) =
+                          Complex (Bx[i], Bz[i]);
+                      }
+
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+            }
+        }
+      else if (typ != MatrixType::Tridiagonal_Hermitian)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+Matrix
+SparseMatrix::bsolve (MatrixType &mattype, const Matrix& b,
+                      octave_idx_type& err, double& rcond,
+                      solve_singularity_handler sing_handler,
+                      bool calc_cond) const
+{
+  Matrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = Matrix (nc, b.cols (), 0.0);
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      volatile int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Banded_Hermitian)
+        {
+          octave_idx_type n_lower = mattype.nlower ();
+          octave_idx_type ldm = n_lower + 1;
+          Matrix m_band (ldm, nc);
+          double *tmp_data = m_band.fortran_vec ();
+
+          if (! mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < ldm; j++)
+                for (octave_idx_type i = 0; i < nc; i++)
+                  tmp_data[ii++] = 0.;
+            }
+
+          for (octave_idx_type j = 0; j < nc; j++)
+            for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+              {
+                octave_idx_type ri = ridx (i);
+                if (ri >= j)
+                  m_band(ri - j, j) = data (i);
+              }
+
+          // Calculate the norm of the matrix, for later use.
+          double anorm;
+          if (calc_cond)
+            anorm = m_band.abs ().sum ().row (0).max ();
+
+          char job = 'L';
+          F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                     nr, n_lower, tmp_data, ldm, err
+                                     F77_CHAR_ARG_LEN (1)));
+
+          if (err != 0)
+            {
+              // Matrix is not positive definite!! Fall through to
+              // unsymmetric banded solver.
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Banded;
+              rcond = 0.0;
+              err = 0;
+            }
+          else
+            {
+              if (calc_cond)
+                {
+                  Array<double> z (dim_vector (3 * nr, 1));
+                  double *pz = z.fortran_vec ();
+                  Array<octave_idx_type> iz (dim_vector (nr, 1));
+                  octave_idx_type *piz = iz.fortran_vec ();
+
+                  F77_XFCN (dpbcon, DPBCON,
+                    (F77_CONST_CHAR_ARG2 (&job, 1),
+                     nr, n_lower, tmp_data, ldm,
+                     anorm, rcond, pz, piz, err
+                     F77_CHAR_ARG_LEN (1)));
+
+                  if (err != 0)
+                    err = -2;
+
+                  volatile double rcond_plus_one = rcond + 1.0;
+
+                  if (rcond_plus_one == 1.0 || xisnan (rcond))
+                    {
+                      err = -2;
+
+                      if (sing_handler)
+                        {
+                          sing_handler (rcond);
+                          mattype.mark_as_rectangular ();
+                        }
+                      else
+                        (*current_liboctave_error_handler)
+                          ("matrix singular to machine precision, rcond = %g",
+                           rcond);
+                    }
+                }
+              else
+                rcond = 1.;
+
+              if (err == 0)
+                {
+                  retval = b;
+                  double *result = retval.fortran_vec ();
+
+                  octave_idx_type b_nc = b.cols ();
+
+                  F77_XFCN (dpbtrs, DPBTRS,
+                            (F77_CONST_CHAR_ARG2 (&job, 1),
+                             nr, n_lower, b_nc, tmp_data,
+                             ldm, result, b.rows (), err
+                             F77_CHAR_ARG_LEN (1)));
+
+                  if (err != 0)
+                    {
+                      (*current_liboctave_error_handler)
+                        ("SparseMatrix::solve solve failed");
+                      err = -1;
+                    }
+                }
+            }
+        }
+
+      if (typ == MatrixType::Banded)
+        {
+          // Create the storage for the banded form of the sparse matrix
+          octave_idx_type n_upper = mattype.nupper ();
+          octave_idx_type n_lower = mattype.nlower ();
+          octave_idx_type ldm = n_upper + 2 * n_lower + 1;
+
+          Matrix m_band (ldm, nc);
+          double *tmp_data = m_band.fortran_vec ();
+
+          if (! mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < ldm; j++)
+                for (octave_idx_type i = 0; i < nc; i++)
+                  tmp_data[ii++] = 0.;
+            }
+
+          for (octave_idx_type j = 0; j < nc; j++)
+            for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+              m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
+
+          // Calculate the norm of the matrix, for later use.
+          double anorm;
+          if (calc_cond)
+            {
+              for (octave_idx_type j = 0; j < nr; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          Array<octave_idx_type> ipvt (dim_vector (nr, 1));
+          octave_idx_type *pipvt = ipvt.fortran_vec ();
+
+          F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
+                                     ldm, pipvt, err));
+
+          // Throw-away extra info LAPACK gives so as to not
+          // change output.
+          if (err != 0)
+            {
+              err = -2;
+              rcond = 0.0;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision");
+
+            }
+          else
+            {
+              if (calc_cond)
+                {
+                  char job = '1';
+                  Array<double> z (dim_vector (3 * nr, 1));
+                  double *pz = z.fortran_vec ();
+                  Array<octave_idx_type> iz (dim_vector (nr, 1));
+                  octave_idx_type *piz = iz.fortran_vec ();
+
+                  F77_XFCN (dgbcon, DGBCON,
+                    (F77_CONST_CHAR_ARG2 (&job, 1),
+                     nc, n_lower, n_upper, tmp_data, ldm, pipvt,
+                     anorm, rcond, pz, piz, err
+                     F77_CHAR_ARG_LEN (1)));
+
+                   if (err != 0)
+                    err = -2;
+
+                  volatile double rcond_plus_one = rcond + 1.0;
+
+                  if (rcond_plus_one == 1.0 || xisnan (rcond))
+                    {
+                      err = -2;
+
+                      if (sing_handler)
+                        {
+                          sing_handler (rcond);
+                          mattype.mark_as_rectangular ();
+                        }
+                      else
+                        (*current_liboctave_error_handler)
+                          ("matrix singular to machine precision, rcond = %g",
+                           rcond);
+                    }
+                }
+              else
+                rcond = 1.;
+
+              if (err == 0)
+                {
+                  retval = b;
+                  double *result = retval.fortran_vec ();
+
+                  octave_idx_type b_nc = b.cols ();
+
+                  char job = 'N';
+                  F77_XFCN (dgbtrs, DGBTRS,
+                            (F77_CONST_CHAR_ARG2 (&job, 1),
+                             nr, n_lower, n_upper, b_nc, tmp_data,
+                             ldm, pipvt, result, b.rows (), err
+                             F77_CHAR_ARG_LEN (1)));
+                }
+            }
+        }
+      else if (typ != MatrixType::Banded_Hermitian)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseMatrix
+SparseMatrix::bsolve (MatrixType &mattype, const SparseMatrix& b,
+                      octave_idx_type& err, double& rcond,
+                      solve_singularity_handler sing_handler,
+                      bool calc_cond) const
+{
+  SparseMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = SparseMatrix (nc, b.cols ());
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      volatile int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Banded_Hermitian)
+        {
+          octave_idx_type n_lower = mattype.nlower ();
+          octave_idx_type ldm = n_lower + 1;
+
+          Matrix m_band (ldm, nc);
+          double *tmp_data = m_band.fortran_vec ();
+
+          if (! mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < ldm; j++)
+                for (octave_idx_type i = 0; i < nc; i++)
+                  tmp_data[ii++] = 0.;
+            }
+
+          for (octave_idx_type j = 0; j < nc; j++)
+            for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+              {
+                octave_idx_type ri = ridx (i);
+                if (ri >= j)
+                  m_band(ri - j, j) = data (i);
+              }
+
+          // Calculate the norm of the matrix, for later use.
+          double anorm;
+          if (calc_cond)
+            anorm = m_band.abs ().sum ().row (0).max ();
+
+          char job = 'L';
+          F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                     nr, n_lower, tmp_data, ldm, err
+                                     F77_CHAR_ARG_LEN (1)));
+
+          if (err != 0)
+            {
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Banded;
+              rcond = 0.0;
+              err = 0;
+            }
+          else
+            {
+              if (calc_cond)
+                {
+                  Array<double> z (dim_vector (3 * nr, 1));
+                  double *pz = z.fortran_vec ();
+                  Array<octave_idx_type> iz (dim_vector (nr, 1));
+                  octave_idx_type *piz = iz.fortran_vec ();
+
+                  F77_XFCN (dpbcon, DPBCON,
+                    (F77_CONST_CHAR_ARG2 (&job, 1),
+                     nr, n_lower, tmp_data, ldm,
+                     anorm, rcond, pz, piz, err
+                     F77_CHAR_ARG_LEN (1)));
+
+                  if (err != 0)
+                    err = -2;
+
+                  volatile double rcond_plus_one = rcond + 1.0;
+
+                  if (rcond_plus_one == 1.0 || xisnan (rcond))
+                    {
+                      err = -2;
+
+                      if (sing_handler)
+                        {
+                          sing_handler (rcond);
+                          mattype.mark_as_rectangular ();
+                        }
+                      else
+                        (*current_liboctave_error_handler)
+                          ("matrix singular to machine precision, rcond = %g",
+                           rcond);
+                    }
+                }
+              else
+                rcond = 1.;
+
+              if (err == 0)
+                {
+                  octave_idx_type b_nr = b.rows ();
+                  octave_idx_type b_nc = b.cols ();
+                  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
+
+                  // Take a first guess that the number of non-zero terms
+                  // will be as many as in b
+                  volatile octave_idx_type x_nz = b.nnz ();
+                  volatile octave_idx_type ii = 0;
+                  retval = SparseMatrix (b_nr, b_nc, x_nz);
+
+                  retval.xcidx (0) = 0;
+                  for (volatile octave_idx_type j = 0; j < b_nc; j++)
+                    {
+                      for (octave_idx_type i = 0; i < b_nr; i++)
+                        Bx[i] = b.elem (i, j);
+
+                      F77_XFCN (dpbtrs, DPBTRS,
+                                (F77_CONST_CHAR_ARG2 (&job, 1),
+                                 nr, n_lower, 1, tmp_data,
+                                 ldm, Bx, b_nr, err
+                                 F77_CHAR_ARG_LEN (1)));
+
+                      if (err != 0)
+                        {
+                          (*current_liboctave_error_handler)
+                            ("SparseMatrix::solve solve failed");
+                          err = -1;
+                          break;
+                        }
+
+                      for (octave_idx_type i = 0; i < b_nr; i++)
+                        {
+                          double tmp = Bx[i];
+                          if (tmp != 0.0)
+                            {
+                              if (ii == x_nz)
+                                {
+                                  // Resize the sparse matrix
+                                  octave_idx_type sz = x_nz *
+                                    (b_nc - j) / b_nc;
+                                  sz = (sz > 10 ? sz : 10) + x_nz;
+                                  retval.change_capacity (sz);
+                                  x_nz = sz;
+                                }
+                              retval.xdata (ii) = tmp;
+                              retval.xridx (ii++) = i;
+                            }
+                        }
+                      retval.xcidx (j+1) = ii;
+                    }
+
+                  retval.maybe_compress ();
+                }
+            }
+        }
+
+      if (typ == MatrixType::Banded)
+        {
+          // Create the storage for the banded form of the sparse matrix
+          octave_idx_type n_upper = mattype.nupper ();
+          octave_idx_type n_lower = mattype.nlower ();
+          octave_idx_type ldm = n_upper + 2 * n_lower + 1;
+
+          Matrix m_band (ldm, nc);
+          double *tmp_data = m_band.fortran_vec ();
+
+          if (! mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < ldm; j++)
+                for (octave_idx_type i = 0; i < nc; i++)
+                  tmp_data[ii++] = 0.;
+            }
+
+          for (octave_idx_type j = 0; j < nc; j++)
+            for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+              m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
+
+          // Calculate the norm of the matrix, for later use.
+          double anorm;
+          if (calc_cond)
+            {
+              for (octave_idx_type j = 0; j < nr; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          Array<octave_idx_type> ipvt (dim_vector (nr, 1));
+          octave_idx_type *pipvt = ipvt.fortran_vec ();
+
+          F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
+                                     ldm, pipvt, err));
+
+          if (err != 0)
+            {
+              err = -2;
+              rcond = 0.0;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision");
+
+            }
+          else
+            {
+              if (calc_cond)
+                {
+                  char job = '1';
+                  Array<double> z (dim_vector (3 * nr, 1));
+                  double *pz = z.fortran_vec ();
+                  Array<octave_idx_type> iz (dim_vector (nr, 1));
+                  octave_idx_type *piz = iz.fortran_vec ();
+
+                  F77_XFCN (dgbcon, DGBCON,
+                    (F77_CONST_CHAR_ARG2 (&job, 1),
+                     nc, n_lower, n_upper, tmp_data, ldm, pipvt,
+                     anorm, rcond, pz, piz, err
+                     F77_CHAR_ARG_LEN (1)));
+
+                   if (err != 0)
+                    err = -2;
+
+                  volatile double rcond_plus_one = rcond + 1.0;
+
+                  if (rcond_plus_one == 1.0 || xisnan (rcond))
+                    {
+                      err = -2;
+
+                      if (sing_handler)
+                        {
+                          sing_handler (rcond);
+                          mattype.mark_as_rectangular ();
+                        }
+                      else
+                        (*current_liboctave_error_handler)
+                          ("matrix singular to machine precision, rcond = %g",
+                           rcond);
+                    }
+                }
+              else
+                rcond = 1.;
+
+              if (err == 0)
+                {
+                  char job = 'N';
+                  volatile octave_idx_type x_nz = b.nnz ();
+                  octave_idx_type b_nc = b.cols ();
+                  retval = SparseMatrix (nr, b_nc, x_nz);
+                  retval.xcidx (0) = 0;
+                  volatile octave_idx_type ii = 0;
+
+                  OCTAVE_LOCAL_BUFFER (double, work, nr);
+
+                  for (volatile octave_idx_type j = 0; j < b_nc; j++)
+                    {
+                      for (octave_idx_type i = 0; i < nr; i++)
+                        work[i] = 0.;
+                      for (octave_idx_type i = b.cidx (j);
+                           i < b.cidx (j+1); i++)
+                        work[b.ridx (i)] = b.data (i);
+
+                      F77_XFCN (dgbtrs, DGBTRS,
+                                (F77_CONST_CHAR_ARG2 (&job, 1),
+                                 nr, n_lower, n_upper, 1, tmp_data,
+                                 ldm, pipvt, work, b.rows (), err
+                                 F77_CHAR_ARG_LEN (1)));
+
+                      // Count non-zeros in work vector and adjust
+                      // space in retval if needed
+                      octave_idx_type new_nnz = 0;
+                      for (octave_idx_type i = 0; i < nr; i++)
+                        if (work[i] != 0.)
+                          new_nnz++;
+
+                      if (ii + new_nnz > x_nz)
+                        {
+                          // Resize the sparse matrix
+                          octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                          retval.change_capacity (sz);
+                          x_nz = sz;
+                        }
+
+                      for (octave_idx_type i = 0; i < nr; i++)
+                        if (work[i] != 0.)
+                          {
+                            retval.xridx (ii) = i;
+                            retval.xdata (ii++) = work[i];
+                          }
+                      retval.xcidx (j+1) = ii;
+                    }
+
+                  retval.maybe_compress ();
+                }
+            }
+        }
+      else if (typ != MatrixType::Banded_Hermitian)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+ComplexMatrix
+SparseMatrix::bsolve (MatrixType &mattype, const ComplexMatrix& b,
+                      octave_idx_type& err, double& rcond,
+                      solve_singularity_handler sing_handler,
+                      bool calc_cond) const
+{
+  ComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      volatile int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Banded_Hermitian)
+        {
+          octave_idx_type n_lower = mattype.nlower ();
+          octave_idx_type ldm = n_lower + 1;
+
+          Matrix m_band (ldm, nc);
+          double *tmp_data = m_band.fortran_vec ();
+
+          if (! mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < ldm; j++)
+                for (octave_idx_type i = 0; i < nc; i++)
+                  tmp_data[ii++] = 0.;
+            }
+
+          for (octave_idx_type j = 0; j < nc; j++)
+            for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+              {
+                octave_idx_type ri = ridx (i);
+                if (ri >= j)
+                  m_band(ri - j, j) = data (i);
+              }
+
+          // Calculate the norm of the matrix, for later use.
+          double anorm;
+          if (calc_cond)
+            anorm = m_band.abs ().sum ().row (0).max ();
+
+          char job = 'L';
+          F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                     nr, n_lower, tmp_data, ldm, err
+                                     F77_CHAR_ARG_LEN (1)));
+
+          if (err != 0)
+            {
+              // Matrix is not positive definite!! Fall through to
+              // unsymmetric banded solver.
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Banded;
+              rcond = 0.0;
+              err = 0;
+            }
+          else
+            {
+              if (calc_cond)
+                {
+                  Array<double> z (dim_vector (3 * nr, 1));
+                  double *pz = z.fortran_vec ();
+                  Array<octave_idx_type> iz (dim_vector (nr, 1));
+                  octave_idx_type *piz = iz.fortran_vec ();
+
+                  F77_XFCN (dpbcon, DPBCON,
+                    (F77_CONST_CHAR_ARG2 (&job, 1),
+                     nr, n_lower, tmp_data, ldm,
+                     anorm, rcond, pz, piz, err
+                     F77_CHAR_ARG_LEN (1)));
+
+                  if (err != 0)
+                    err = -2;
+
+                  volatile double rcond_plus_one = rcond + 1.0;
+
+                  if (rcond_plus_one == 1.0 || xisnan (rcond))
+                    {
+                      err = -2;
+
+                      if (sing_handler)
+                        {
+                          sing_handler (rcond);
+                          mattype.mark_as_rectangular ();
+                        }
+                      else
+                        (*current_liboctave_error_handler)
+                          ("matrix singular to machine precision, rcond = %g",
+                           rcond);
+                    }
+                }
+              else
+                rcond = 1.;
+
+              if (err == 0)
+                {
+                  octave_idx_type b_nr = b.rows ();
+                  octave_idx_type b_nc = b.cols ();
+
+                  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
+                  OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
+
+                  retval.resize (b_nr, b_nc);
+
+                  for (volatile octave_idx_type j = 0; j < b_nc; j++)
+                    {
+                      for (octave_idx_type i = 0; i < b_nr; i++)
+                        {
+                          Complex c = b (i,j);
+                          Bx[i] = std::real (c);
+                          Bz[i] = std::imag (c);
+                        }
+
+                      F77_XFCN (dpbtrs, DPBTRS,
+                                (F77_CONST_CHAR_ARG2 (&job, 1),
+                                 nr, n_lower, 1, tmp_data,
+                                 ldm, Bx, b_nr, err
+                                 F77_CHAR_ARG_LEN (1)));
+
+                      if (err != 0)
+                        {
+                          (*current_liboctave_error_handler)
+                            ("SparseMatrix::solve solve failed");
+                          err = -1;
+                          break;
+                        }
+
+                      F77_XFCN (dpbtrs, DPBTRS,
+                                (F77_CONST_CHAR_ARG2 (&job, 1),
+                                 nr, n_lower, 1, tmp_data,
+                                 ldm, Bz, b.rows (), err
+                                 F77_CHAR_ARG_LEN (1)));
+
+                      if (err != 0)
+                        {
+                          (*current_liboctave_error_handler)
+                            ("SparseMatrix::solve solve failed");
+                          err = -1;
+                          break;
+                        }
+
+                      for (octave_idx_type i = 0; i < b_nr; i++)
+                        retval(i, j) = Complex (Bx[i], Bz[i]);
+                    }
+                }
+            }
+        }
+
+      if (typ == MatrixType::Banded)
+        {
+          // Create the storage for the banded form of the sparse matrix
+          octave_idx_type n_upper = mattype.nupper ();
+          octave_idx_type n_lower = mattype.nlower ();
+          octave_idx_type ldm = n_upper + 2 * n_lower + 1;
+
+          Matrix m_band (ldm, nc);
+          double *tmp_data = m_band.fortran_vec ();
+
+          if (! mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < ldm; j++)
+                for (octave_idx_type i = 0; i < nc; i++)
+                  tmp_data[ii++] = 0.;
+            }
+
+          for (octave_idx_type j = 0; j < nc; j++)
+            for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+              m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
+
+          // Calculate the norm of the matrix, for later use.
+          double anorm;
+          if (calc_cond)
+            {
+              for (octave_idx_type j = 0; j < nr; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          Array<octave_idx_type> ipvt (dim_vector (nr, 1));
+          octave_idx_type *pipvt = ipvt.fortran_vec ();
+
+          F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
+                                     ldm, pipvt, err));
+
+          if (err != 0)
+            {
+              err = -2;
+              rcond = 0.0;
+
+              if (sing_handler)
+                {
+                sing_handler (rcond);
+                mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision");
+
+            }
+          else
+            {
+              if (calc_cond)
+                {
+                  char job = '1';
+                  Array<double> z (dim_vector (3 * nr, 1));
+                  double *pz = z.fortran_vec ();
+                  Array<octave_idx_type> iz (dim_vector (nr, 1));
+                  octave_idx_type *piz = iz.fortran_vec ();
+
+                  F77_XFCN (dpbcon, DPBCON,
+                    (F77_CONST_CHAR_ARG2 (&job, 1),
+                     nr, n_lower, tmp_data, ldm,
+                     anorm, rcond, pz, piz, err
+                     F77_CHAR_ARG_LEN (1)));
+
+                  if (err != 0)
+                    err = -2;
+
+                  volatile double rcond_plus_one = rcond + 1.0;
+
+                  if (rcond_plus_one == 1.0 || xisnan (rcond))
+                    {
+                      err = -2;
+
+                      if (sing_handler)
+                        {
+                        sing_handler (rcond);
+                        mattype.mark_as_rectangular ();
+                        }
+                      else
+                        (*current_liboctave_error_handler)
+                          ("matrix singular to machine precision, rcond = %g",
+                           rcond);
+                    }
+                }
+              else
+                rcond = 1.;
+
+              if (err == 0)
+                {
+                  char job = 'N';
+                  octave_idx_type b_nc = b.cols ();
+                  retval.resize (nr,b_nc);
+
+                  OCTAVE_LOCAL_BUFFER (double, Bz, nr);
+                  OCTAVE_LOCAL_BUFFER (double, Bx, nr);
+
+                  for (volatile octave_idx_type j = 0; j < b_nc; j++)
+                    {
+                      for (octave_idx_type i = 0; i < nr; i++)
+                        {
+                          Complex c = b (i, j);
+                          Bx[i] = std::real (c);
+                          Bz[i] = std::imag  (c);
+                        }
+
+                      F77_XFCN (dgbtrs, DGBTRS,
+                                (F77_CONST_CHAR_ARG2 (&job, 1),
+                                 nr, n_lower, n_upper, 1, tmp_data,
+                                 ldm, pipvt, Bx, b.rows (), err
+                                 F77_CHAR_ARG_LEN (1)));
+
+                      F77_XFCN (dgbtrs, DGBTRS,
+                                (F77_CONST_CHAR_ARG2 (&job, 1),
+                                 nr, n_lower, n_upper, 1, tmp_data,
+                                 ldm, pipvt, Bz, b.rows (), err
+                                 F77_CHAR_ARG_LEN (1)));
+
+                      for (octave_idx_type i = 0; i < nr; i++)
+                        retval(i, j) = Complex (Bx[i], Bz[i]);
+                    }
+                }
+            }
+        }
+      else if (typ != MatrixType::Banded_Hermitian)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseComplexMatrix
+SparseMatrix::bsolve (MatrixType &mattype, const SparseComplexMatrix& b,
+                      octave_idx_type& err, double& rcond,
+                      solve_singularity_handler sing_handler,
+                      bool calc_cond) const
+{
+  SparseComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = SparseComplexMatrix (nc, b.cols ());
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      volatile int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Banded_Hermitian)
+        {
+          octave_idx_type n_lower = mattype.nlower ();
+          octave_idx_type ldm = n_lower + 1;
+
+          Matrix m_band (ldm, nc);
+          double *tmp_data = m_band.fortran_vec ();
+
+          if (! mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < ldm; j++)
+                for (octave_idx_type i = 0; i < nc; i++)
+                  tmp_data[ii++] = 0.;
+            }
+
+          for (octave_idx_type j = 0; j < nc; j++)
+            for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+              {
+                octave_idx_type ri = ridx (i);
+                if (ri >= j)
+                  m_band(ri - j, j) = data (i);
+              }
+
+          // Calculate the norm of the matrix, for later use.
+          double anorm;
+          if (calc_cond)
+            anorm = m_band.abs ().sum ().row (0).max ();
+
+          char job = 'L';
+          F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
+                                     nr, n_lower, tmp_data, ldm, err
+                                     F77_CHAR_ARG_LEN (1)));
+
+          if (err != 0)
+            {
+              // Matrix is not positive definite!! Fall through to
+              // unsymmetric banded solver.
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Banded;
+
+              rcond = 0.0;
+              err = 0;
+            }
+          else
+            {
+              if (calc_cond)
+                {
+                  Array<double> z (dim_vector (3 * nr, 1));
+                  double *pz = z.fortran_vec ();
+                  Array<octave_idx_type> iz (dim_vector (nr, 1));
+                  octave_idx_type *piz = iz.fortran_vec ();
+
+                  F77_XFCN (dpbcon, DPBCON,
+                    (F77_CONST_CHAR_ARG2 (&job, 1),
+                     nr, n_lower, tmp_data, ldm,
+                     anorm, rcond, pz, piz, err
+                     F77_CHAR_ARG_LEN (1)));
+
+                  if (err != 0)
+                    err = -2;
+
+                  volatile double rcond_plus_one = rcond + 1.0;
+
+                  if (rcond_plus_one == 1.0 || xisnan (rcond))
+                    {
+                      err = -2;
+
+                      if (sing_handler)
+                        {
+                          sing_handler (rcond);
+                          mattype.mark_as_rectangular ();
+                        }
+                      else
+                        (*current_liboctave_error_handler)
+                          ("matrix singular to machine precision, rcond = %g",
+                           rcond);
+                    }
+                }
+              else
+                rcond = 1.;
+
+              if (err == 0)
+                {
+                  octave_idx_type b_nr = b.rows ();
+                  octave_idx_type b_nc = b.cols ();
+                  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
+                  OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
+
+                  // Take a first guess that the number of non-zero terms
+                  // will be as many as in b
+                  volatile octave_idx_type x_nz = b.nnz ();
+                  volatile octave_idx_type ii = 0;
+                  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
+
+                  retval.xcidx (0) = 0;
+                  for (volatile octave_idx_type j = 0; j < b_nc; j++)
+                    {
+
+                      for (octave_idx_type i = 0; i < b_nr; i++)
+                        {
+                          Complex c = b (i,j);
+                          Bx[i] = std::real (c);
+                          Bz[i] = std::imag (c);
+                        }
+
+                      F77_XFCN (dpbtrs, DPBTRS,
+                                (F77_CONST_CHAR_ARG2 (&job, 1),
+                                 nr, n_lower, 1, tmp_data,
+                                 ldm, Bx, b_nr, err
+                                 F77_CHAR_ARG_LEN (1)));
+
+                      if (err != 0)
+                        {
+                          (*current_liboctave_error_handler)
+                            ("SparseMatrix::solve solve failed");
+                          err = -1;
+                          break;
+                        }
+
+                      F77_XFCN (dpbtrs, DPBTRS,
+                                (F77_CONST_CHAR_ARG2 (&job, 1),
+                                 nr, n_lower, 1, tmp_data,
+                                 ldm, Bz, b_nr, err
+                                 F77_CHAR_ARG_LEN (1)));
+
+                      if (err != 0)
+                        {
+                          (*current_liboctave_error_handler)
+                            ("SparseMatrix::solve solve failed");
+
+                          err = -1;
+                          break;
+                        }
+
+                      // Count non-zeros in work vector and adjust
+                      // space in retval if needed
+                      octave_idx_type new_nnz = 0;
+                      for (octave_idx_type i = 0; i < nr; i++)
+                        if (Bx[i] != 0. || Bz[i] != 0.)
+                          new_nnz++;
+
+                      if (ii + new_nnz > x_nz)
+                        {
+                          // Resize the sparse matrix
+                          octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                          retval.change_capacity (sz);
+                          x_nz = sz;
+                        }
+
+                      for (octave_idx_type i = 0; i < nr; i++)
+                        if (Bx[i] != 0. || Bz[i] != 0.)
+                          {
+                            retval.xridx (ii) = i;
+                            retval.xdata (ii++) =
+                              Complex (Bx[i], Bz[i]);
+                          }
+
+                      retval.xcidx (j+1) = ii;
+                    }
+
+                  retval.maybe_compress ();
+                }
+            }
+        }
+
+      if (typ == MatrixType::Banded)
+        {
+          // Create the storage for the banded form of the sparse matrix
+          octave_idx_type n_upper = mattype.nupper ();
+          octave_idx_type n_lower = mattype.nlower ();
+          octave_idx_type ldm = n_upper + 2 * n_lower + 1;
+
+          Matrix m_band (ldm, nc);
+          double *tmp_data = m_band.fortran_vec ();
+
+          if (! mattype.is_dense ())
+            {
+              octave_idx_type ii = 0;
+
+              for (octave_idx_type j = 0; j < ldm; j++)
+                for (octave_idx_type i = 0; i < nc; i++)
+                  tmp_data[ii++] = 0.;
+            }
+
+          for (octave_idx_type j = 0; j < nc; j++)
+            for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+              m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
+
+          // Calculate the norm of the matrix, for later use.
+          double anorm;
+          if (calc_cond)
+            {
+              for (octave_idx_type j = 0; j < nr; j++)
+                {
+                  double atmp = 0.;
+                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                    atmp += fabs (data (i));
+                  if (atmp > anorm)
+                    anorm = atmp;
+                }
+            }
+
+          Array<octave_idx_type> ipvt (dim_vector (nr, 1));
+          octave_idx_type *pipvt = ipvt.fortran_vec ();
+
+          F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
+                                     ldm, pipvt, err));
+
+          if (err != 0)
+            {
+              err = -2;
+              rcond = 0.0;
+
+              if (sing_handler)
+                {
+                  sing_handler (rcond);
+                  mattype.mark_as_rectangular ();
+                }
+              else
+                (*current_liboctave_error_handler)
+                  ("matrix singular to machine precision");
+
+            }
+          else
+            {
+              if (calc_cond)
+                {
+                  char job = '1';
+                  Array<double> z (dim_vector (3 * nr, 1));
+                  double *pz = z.fortran_vec ();
+                  Array<octave_idx_type> iz (dim_vector (nr, 1));
+                  octave_idx_type *piz = iz.fortran_vec ();
+
+                  F77_XFCN (dgbcon, DGBCON,
+                    (F77_CONST_CHAR_ARG2 (&job, 1),
+                     nc, n_lower, n_upper, tmp_data, ldm, pipvt,
+                     anorm, rcond, pz, piz, err
+                     F77_CHAR_ARG_LEN (1)));
+
+                   if (err != 0)
+                    err = -2;
+
+                  volatile double rcond_plus_one = rcond + 1.0;
+
+                  if (rcond_plus_one == 1.0 || xisnan (rcond))
+                    {
+                      err = -2;
+
+                      if (sing_handler)
+                        {
+                          sing_handler (rcond);
+                          mattype.mark_as_rectangular ();
+                        }
+                      else
+                        (*current_liboctave_error_handler)
+                          ("matrix singular to machine precision, rcond = %g",
+                           rcond);
+                    }
+                }
+              else
+                rcond = 1.;
+
+              if (err == 0)
+                {
+                  char job = 'N';
+                  volatile octave_idx_type x_nz = b.nnz ();
+                  octave_idx_type b_nc = b.cols ();
+                  retval = SparseComplexMatrix (nr, b_nc, x_nz);
+                  retval.xcidx (0) = 0;
+                  volatile octave_idx_type ii = 0;
+
+                  OCTAVE_LOCAL_BUFFER (double, Bx, nr);
+                  OCTAVE_LOCAL_BUFFER (double, Bz, nr);
+
+                  for (volatile octave_idx_type j = 0; j < b_nc; j++)
+                    {
+                      for (octave_idx_type i = 0; i < nr; i++)
+                        {
+                          Bx[i] = 0.;
+                          Bz[i] = 0.;
+                        }
+                      for (octave_idx_type i = b.cidx (j);
+                           i < b.cidx (j+1); i++)
+                        {
+                          Complex c = b.data (i);
+                          Bx[b.ridx (i)] = std::real (c);
+                          Bz[b.ridx (i)] = std::imag (c);
+                        }
+
+                      F77_XFCN (dgbtrs, DGBTRS,
+                                (F77_CONST_CHAR_ARG2 (&job, 1),
+                                 nr, n_lower, n_upper, 1, tmp_data,
+                                 ldm, pipvt, Bx, b.rows (), err
+                                 F77_CHAR_ARG_LEN (1)));
+
+                      F77_XFCN (dgbtrs, DGBTRS,
+                                (F77_CONST_CHAR_ARG2 (&job, 1),
+                                 nr, n_lower, n_upper, 1, tmp_data,
+                                 ldm, pipvt, Bz, b.rows (), err
+                                 F77_CHAR_ARG_LEN (1)));
+
+                      // Count non-zeros in work vector and adjust
+                      // space in retval if needed
+                      octave_idx_type new_nnz = 0;
+                      for (octave_idx_type i = 0; i < nr; i++)
+                        if (Bx[i] != 0. || Bz[i] != 0.)
+                          new_nnz++;
+
+                      if (ii + new_nnz > x_nz)
+                        {
+                          // Resize the sparse matrix
+                          octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                          retval.change_capacity (sz);
+                          x_nz = sz;
+                        }
+
+                      for (octave_idx_type i = 0; i < nr; i++)
+                        if (Bx[i] != 0. || Bz[i] != 0.)
+                          {
+                            retval.xridx (ii) = i;
+                            retval.xdata (ii++) =
+                              Complex (Bx[i], Bz[i]);
+                          }
+                      retval.xcidx (j+1) = ii;
+                    }
+
+                  retval.maybe_compress ();
+                }
+            }
+        }
+      else if (typ != MatrixType::Banded_Hermitian)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+void *
+SparseMatrix::factorize (octave_idx_type& err, double &rcond, Matrix &Control,
+                         Matrix &Info, solve_singularity_handler sing_handler,
+                         bool calc_cond) const
+{
+  // The return values
+  void *Numeric = 0;
+  err = 0;
+
+#ifdef HAVE_UMFPACK
+  // Setup the control parameters
+  Control = Matrix (UMFPACK_CONTROL, 1);
+  double *control = Control.fortran_vec ();
+  UMFPACK_DNAME (defaults) (control);
+
+  double tmp = octave_sparse_params::get_key ("spumoni");
+  if (!xisnan (tmp))
+    Control (UMFPACK_PRL) = tmp;
+  tmp = octave_sparse_params::get_key ("piv_tol");
+  if (!xisnan (tmp))
+    {
+      Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
+      Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
+    }
+
+  // Set whether we are allowed to modify Q or not
+  tmp = octave_sparse_params::get_key ("autoamd");
+  if (!xisnan (tmp))
+    Control (UMFPACK_FIXQ) = tmp;
+
+  UMFPACK_DNAME (report_control) (control);
+
+  const octave_idx_type *Ap = cidx ();
+  const octave_idx_type *Ai = ridx ();
+  const double *Ax = data ();
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+
+  UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1, control);
+
+  void *Symbolic;
+  Info = Matrix (1, UMFPACK_INFO);
+  double *info = Info.fortran_vec ();
+  int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, 0,
+                                     &Symbolic, control, info);
+
+  if (status < 0)
+    {
+      (*current_liboctave_error_handler)
+        ("SparseMatrix::solve symbolic factorization failed");
+      err = -1;
+
+      UMFPACK_DNAME (report_status) (control, status);
+      UMFPACK_DNAME (report_info) (control, info);
+
+      UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
+    }
+  else
+    {
+      UMFPACK_DNAME (report_symbolic) (Symbolic, control);
+
+      status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic,
+                                   &Numeric, control, info) ;
+      UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
+
+      if (calc_cond)
+        rcond = Info (UMFPACK_RCOND);
+      else
+        rcond = 1.;
+      volatile double rcond_plus_one = rcond + 1.0;
+
+      if (status == UMFPACK_WARNING_singular_matrix ||
+          rcond_plus_one == 1.0 || xisnan (rcond))
+        {
+          UMFPACK_DNAME (report_numeric) (Numeric, control);
+
+          err = -2;
+
+          if (sing_handler)
+            sing_handler (rcond);
+          else
+            (*current_liboctave_error_handler)
+              ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+               rcond);
+
+        }
+      else if (status < 0)
+          {
+            (*current_liboctave_error_handler)
+              ("SparseMatrix::solve numeric factorization failed");
+
+            UMFPACK_DNAME (report_status) (control, status);
+            UMFPACK_DNAME (report_info) (control, info);
+
+            err = -1;
+          }
+        else
+          {
+            UMFPACK_DNAME (report_numeric) (Numeric, control);
+          }
+    }
+
+  if (err != 0)
+    UMFPACK_DNAME (free_numeric) (&Numeric);
+
+#else
+  (*current_liboctave_error_handler) ("UMFPACK not installed");
+#endif
+
+  return Numeric;
+}
+
+Matrix
+SparseMatrix::fsolve (MatrixType &mattype, const Matrix& b,
+                      octave_idx_type& err, double& rcond,
+                      solve_singularity_handler sing_handler,
+                      bool calc_cond) const
+{
+  Matrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = Matrix (nc, b.cols (), 0.0);
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      volatile int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Hermitian)
+        {
+#ifdef HAVE_CHOLMOD
+          cholmod_common Common;
+          cholmod_common *cm = &Common;
+
+          // Setup initial parameters
+          CHOLMOD_NAME(start) (cm);
+          cm->prefer_zomplex = false;
+
+          double spu = octave_sparse_params::get_key ("spumoni");
+          if (spu == 0.)
+            {
+              cm->print = -1;
+              cm->print_function = 0;
+            }
+          else
+            {
+              cm->print = static_cast<int> (spu) + 2;
+              cm->print_function =&SparseCholPrint;
+            }
+
+          cm->error_handler = &SparseCholError;
+          cm->complex_divide = CHOLMOD_NAME(divcomplex);
+          cm->hypotenuse = CHOLMOD_NAME(hypot);
+
+          cm->final_ll = true;
+
+          cholmod_sparse Astore;
+          cholmod_sparse *A = &Astore;
+          double dummy;
+          A->nrow = nr;
+          A->ncol = nc;
+
+          A->p = cidx ();
+          A->i = ridx ();
+          A->nzmax = nnz ();
+          A->packed = true;
+          A->sorted = true;
+          A->nz = 0;
+#ifdef IDX_TYPE_LONG
+          A->itype = CHOLMOD_LONG;
+#else
+          A->itype = CHOLMOD_INT;
+#endif
+          A->dtype = CHOLMOD_DOUBLE;
+          A->stype = 1;
+          A->xtype = CHOLMOD_REAL;
+
+          if (nr < 1)
+            A->x = &dummy;
+          else
+            A->x = data ();
+
+          cholmod_dense Bstore;
+          cholmod_dense *B = &Bstore;
+          B->nrow = b.rows ();
+          B->ncol = b.cols ();
+          B->d = B->nrow;
+          B->nzmax = B->nrow * B->ncol;
+          B->dtype = CHOLMOD_DOUBLE;
+          B->xtype = CHOLMOD_REAL;
+          if (nc < 1 || b.cols () < 1)
+            B->x = &dummy;
+          else
+            // We won't alter it, honest :-)
+            B->x = const_cast<double *>(b.fortran_vec ());
+
+          cholmod_factor *L;
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          L = CHOLMOD_NAME(analyze) (A, cm);
+          CHOLMOD_NAME(factorize) (A, L, cm);
+          if (calc_cond)
+            rcond = CHOLMOD_NAME(rcond)(L, cm);
+          else
+            rcond = 1.0;
+
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+          if (rcond == 0.0)
+            {
+              // Either its indefinite or singular. Try UMFPACK
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Full;
+            }
+          else
+            {
+              volatile double rcond_plus_one = rcond + 1.0;
+
+              if (rcond_plus_one == 1.0 || xisnan (rcond))
+                {
+                  err = -2;
+
+                  if (sing_handler)
+                    {
+                      sing_handler (rcond);
+                      mattype.mark_as_rectangular ();
+                    }
+                  else
+                    (*current_liboctave_error_handler)
+                      ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                       rcond);
+
+                  return retval;
+                }
+
+              cholmod_dense *X;
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+              retval.resize (b.rows (), b.cols ());
+              for (octave_idx_type j = 0; j < b.cols (); j++)
+                {
+                  octave_idx_type jr = j * b.rows ();
+                  for (octave_idx_type i = 0; i < b.rows (); i++)
+                    retval.xelem (i,j) = static_cast<double *>(X->x)[jr + i];
+                }
+
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              CHOLMOD_NAME(free_dense) (&X, cm);
+              CHOLMOD_NAME(free_factor) (&L, cm);
+              CHOLMOD_NAME(finish) (cm);
+              static char tmp[] = " ";
+              CHOLMOD_NAME(print_common) (tmp, cm);
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+            }
+#else
+          (*current_liboctave_warning_handler)
+            ("CHOLMOD not installed");
+
+          mattype.mark_as_unsymmetric ();
+          typ = MatrixType::Full;
+#endif
+        }
+
+      if (typ == MatrixType::Full)
+        {
+#ifdef HAVE_UMFPACK
+          Matrix Control, Info;
+          void *Numeric =
+            factorize (err, rcond, Control, Info, sing_handler, calc_cond);
+
+          if (err == 0)
+            {
+              const double *Bx = b.fortran_vec ();
+              retval.resize (b.rows (), b.cols ());
+              double *result = retval.fortran_vec ();
+              octave_idx_type b_nr = b.rows ();
+              octave_idx_type b_nc = b.cols ();
+              int status = 0;
+              double *control = Control.fortran_vec ();
+              double *info = Info.fortran_vec ();
+              const octave_idx_type *Ap = cidx ();
+              const octave_idx_type *Ai = ridx ();
+              const double *Ax = data ();
+
+              for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
+                {
+                  status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap,
+                                             Ai, Ax, &result[iidx], &Bx[iidx],
+                                             Numeric, control, info);
+                  if (status < 0)
+                    {
+                      (*current_liboctave_error_handler)
+                        ("SparseMatrix::solve solve failed");
+
+                      UMFPACK_DNAME (report_status) (control, status);
+
+                      err = -1;
+
+                      break;
+                    }
+                }
+
+              UMFPACK_DNAME (report_info) (control, info);
+
+              UMFPACK_DNAME (free_numeric) (&Numeric);
+            }
+          else
+            mattype.mark_as_rectangular ();
+
+#else
+          (*current_liboctave_error_handler) ("UMFPACK not installed");
+#endif
+        }
+      else if (typ != MatrixType::Hermitian)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseMatrix
+SparseMatrix::fsolve (MatrixType &mattype, const SparseMatrix& b,
+                      octave_idx_type& err, double& rcond,
+                      solve_singularity_handler sing_handler,
+                      bool calc_cond) const
+{
+  SparseMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = SparseMatrix (nc, b.cols ());
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      volatile int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Hermitian)
+        {
+#ifdef HAVE_CHOLMOD
+          cholmod_common Common;
+          cholmod_common *cm = &Common;
+
+          // Setup initial parameters
+          CHOLMOD_NAME(start) (cm);
+          cm->prefer_zomplex = false;
+
+          double spu = octave_sparse_params::get_key ("spumoni");
+          if (spu == 0.)
+            {
+              cm->print = -1;
+              cm->print_function = 0;
+            }
+          else
+            {
+              cm->print = static_cast<int> (spu) + 2;
+              cm->print_function =&SparseCholPrint;
+            }
+
+          cm->error_handler = &SparseCholError;
+          cm->complex_divide = CHOLMOD_NAME(divcomplex);
+          cm->hypotenuse = CHOLMOD_NAME(hypot);
+
+          cm->final_ll = true;
+
+          cholmod_sparse Astore;
+          cholmod_sparse *A = &Astore;
+          double dummy;
+          A->nrow = nr;
+          A->ncol = nc;
+
+          A->p = cidx ();
+          A->i = ridx ();
+          A->nzmax = nnz ();
+          A->packed = true;
+          A->sorted = true;
+          A->nz = 0;
+#ifdef IDX_TYPE_LONG
+          A->itype = CHOLMOD_LONG;
+#else
+          A->itype = CHOLMOD_INT;
+#endif
+          A->dtype = CHOLMOD_DOUBLE;
+          A->stype = 1;
+          A->xtype = CHOLMOD_REAL;
+
+          if (nr < 1)
+            A->x = &dummy;
+          else
+            A->x = data ();
+
+          cholmod_sparse Bstore;
+          cholmod_sparse *B = &Bstore;
+          B->nrow = b.rows ();
+          B->ncol = b.cols ();
+          B->p = b.cidx ();
+          B->i = b.ridx ();
+          B->nzmax = b.nnz ();
+          B->packed = true;
+          B->sorted = true;
+          B->nz = 0;
+#ifdef IDX_TYPE_LONG
+          B->itype = CHOLMOD_LONG;
+#else
+          B->itype = CHOLMOD_INT;
+#endif
+          B->dtype = CHOLMOD_DOUBLE;
+          B->stype = 0;
+          B->xtype = CHOLMOD_REAL;
+
+          if (b.rows () < 1 || b.cols () < 1)
+            B->x = &dummy;
+          else
+            B->x = b.data ();
+
+          cholmod_factor *L;
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          L = CHOLMOD_NAME(analyze) (A, cm);
+          CHOLMOD_NAME(factorize) (A, L, cm);
+          if (calc_cond)
+            rcond = CHOLMOD_NAME(rcond)(L, cm);
+          else
+            rcond = 1.;
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+          if (rcond == 0.0)
+            {
+              // Either its indefinite or singular. Try UMFPACK
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Full;
+            }
+          else
+            {
+              volatile double rcond_plus_one = rcond + 1.0;
+
+              if (rcond_plus_one == 1.0 || xisnan (rcond))
+                {
+                  err = -2;
+
+                  if (sing_handler)
+                    {
+                      sing_handler (rcond);
+                      mattype.mark_as_rectangular ();
+                    }
+                  else
+                    (*current_liboctave_error_handler)
+                      ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                       rcond);
+
+                  return retval;
+                }
+
+              cholmod_sparse *X;
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+              retval = SparseMatrix (static_cast<octave_idx_type>(X->nrow),
+                                     static_cast<octave_idx_type>(X->ncol),
+                                     static_cast<octave_idx_type>(X->nzmax));
+              for (octave_idx_type j = 0;
+                   j <= static_cast<octave_idx_type>(X->ncol); j++)
+                retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
+              for (octave_idx_type j = 0;
+                   j < static_cast<octave_idx_type>(X->nzmax); j++)
+                {
+                  retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
+                  retval.xdata (j) = static_cast<double *>(X->x)[j];
+                }
+
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              CHOLMOD_NAME(free_sparse) (&X, cm);
+              CHOLMOD_NAME(free_factor) (&L, cm);
+              CHOLMOD_NAME(finish) (cm);
+              static char tmp[] = " ";
+              CHOLMOD_NAME(print_common) (tmp, cm);
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+            }
+#else
+          (*current_liboctave_warning_handler)
+            ("CHOLMOD not installed");
+
+          mattype.mark_as_unsymmetric ();
+          typ = MatrixType::Full;
+#endif
+        }
+
+      if (typ == MatrixType::Full)
+        {
+#ifdef HAVE_UMFPACK
+          Matrix Control, Info;
+          void *Numeric = factorize (err, rcond, Control, Info,
+                                     sing_handler, calc_cond);
+
+          if (err == 0)
+            {
+              octave_idx_type b_nr = b.rows ();
+              octave_idx_type b_nc = b.cols ();
+              int status = 0;
+              double *control = Control.fortran_vec ();
+              double *info = Info.fortran_vec ();
+              const octave_idx_type *Ap = cidx ();
+              const octave_idx_type *Ai = ridx ();
+              const double *Ax = data ();
+
+              OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
+              OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
+
+              // Take a first guess that the number of non-zero terms
+              // will be as many as in b
+              octave_idx_type x_nz = b.nnz ();
+              octave_idx_type ii = 0;
+              retval = SparseMatrix (b_nr, b_nc, x_nz);
+
+              retval.xcidx (0) = 0;
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+
+                  for (octave_idx_type i = 0; i < b_nr; i++)
+                    Bx[i] = b.elem (i, j);
+
+                  status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap,
+                                             Ai, Ax, Xx, Bx, Numeric, control,
+                                             info);
+                  if (status < 0)
+                    {
+                      (*current_liboctave_error_handler)
+                        ("SparseMatrix::solve solve failed");
+
+                      UMFPACK_DNAME (report_status) (control, status);
+
+                      err = -1;
+
+                      break;
+                    }
+
+                  for (octave_idx_type i = 0; i < b_nr; i++)
+                    {
+                      double tmp = Xx[i];
+                      if (tmp != 0.0)
+                        {
+                          if (ii == x_nz)
+                            {
+                              // Resize the sparse matrix
+                              octave_idx_type sz = x_nz * (b_nc - j) / b_nc;
+                              sz = (sz > 10 ? sz : 10) + x_nz;
+                              retval.change_capacity (sz);
+                              x_nz = sz;
+                            }
+                          retval.xdata (ii) = tmp;
+                          retval.xridx (ii++) = i;
+                        }
+                    }
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+
+              UMFPACK_DNAME (report_info) (control, info);
+
+              UMFPACK_DNAME (free_numeric) (&Numeric);
+            }
+          else
+            mattype.mark_as_rectangular ();
+
+#else
+          (*current_liboctave_error_handler) ("UMFPACK not installed");
+#endif
+        }
+      else if (typ != MatrixType::Hermitian)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+ComplexMatrix
+SparseMatrix::fsolve (MatrixType &mattype, const ComplexMatrix& b,
+                      octave_idx_type& err, double& rcond,
+                      solve_singularity_handler sing_handler,
+                      bool calc_cond) const
+{
+  ComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      volatile int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Hermitian)
+        {
+#ifdef HAVE_CHOLMOD
+          cholmod_common Common;
+          cholmod_common *cm = &Common;
+
+          // Setup initial parameters
+          CHOLMOD_NAME(start) (cm);
+          cm->prefer_zomplex = false;
+
+          double spu = octave_sparse_params::get_key ("spumoni");
+          if (spu == 0.)
+            {
+              cm->print = -1;
+              cm->print_function = 0;
+            }
+          else
+            {
+              cm->print = static_cast<int> (spu) + 2;
+              cm->print_function =&SparseCholPrint;
+            }
+
+          cm->error_handler = &SparseCholError;
+          cm->complex_divide = CHOLMOD_NAME(divcomplex);
+          cm->hypotenuse = CHOLMOD_NAME(hypot);
+
+          cm->final_ll = true;
+
+          cholmod_sparse Astore;
+          cholmod_sparse *A = &Astore;
+          double dummy;
+          A->nrow = nr;
+          A->ncol = nc;
+
+          A->p = cidx ();
+          A->i = ridx ();
+          A->nzmax = nnz ();
+          A->packed = true;
+          A->sorted = true;
+          A->nz = 0;
+#ifdef IDX_TYPE_LONG
+          A->itype = CHOLMOD_LONG;
+#else
+          A->itype = CHOLMOD_INT;
+#endif
+          A->dtype = CHOLMOD_DOUBLE;
+          A->stype = 1;
+          A->xtype = CHOLMOD_REAL;
+
+          if (nr < 1)
+            A->x = &dummy;
+          else
+            A->x = data ();
+
+          cholmod_dense Bstore;
+          cholmod_dense *B = &Bstore;
+          B->nrow = b.rows ();
+          B->ncol = b.cols ();
+          B->d = B->nrow;
+          B->nzmax = B->nrow * B->ncol;
+          B->dtype = CHOLMOD_DOUBLE;
+          B->xtype = CHOLMOD_COMPLEX;
+          if (nc < 1 || b.cols () < 1)
+            B->x = &dummy;
+          else
+            // We won't alter it, honest :-)
+            B->x = const_cast<Complex *>(b.fortran_vec ());
+
+          cholmod_factor *L;
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          L = CHOLMOD_NAME(analyze) (A, cm);
+          CHOLMOD_NAME(factorize) (A, L, cm);
+          if (calc_cond)
+            rcond = CHOLMOD_NAME(rcond)(L, cm);
+          else
+            rcond = 1.0;
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+          if (rcond == 0.0)
+            {
+              // Either its indefinite or singular. Try UMFPACK
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Full;
+            }
+          else
+            {
+              volatile double rcond_plus_one = rcond + 1.0;
+
+              if (rcond_plus_one == 1.0 || xisnan (rcond))
+                {
+                  err = -2;
+
+                  if (sing_handler)
+                    {
+                      sing_handler (rcond);
+                      mattype.mark_as_rectangular ();
+                    }
+                  else
+                    (*current_liboctave_error_handler)
+                      ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                       rcond);
+
+                  return retval;
+                }
+
+              cholmod_dense *X;
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+              retval.resize (b.rows (), b.cols ());
+              for (octave_idx_type j = 0; j < b.cols (); j++)
+                {
+                  octave_idx_type jr = j * b.rows ();
+                  for (octave_idx_type i = 0; i < b.rows (); i++)
+                    retval.xelem (i,j) = static_cast<Complex *>(X->x)[jr + i];
+                }
+
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              CHOLMOD_NAME(free_dense) (&X, cm);
+              CHOLMOD_NAME(free_factor) (&L, cm);
+              CHOLMOD_NAME(finish) (cm);
+              static char tmp[] = " ";
+              CHOLMOD_NAME(print_common) (tmp, cm);
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+            }
+#else
+          (*current_liboctave_warning_handler)
+            ("CHOLMOD not installed");
+
+          mattype.mark_as_unsymmetric ();
+          typ = MatrixType::Full;
+#endif
+        }
+
+      if (typ == MatrixType::Full)
+        {
+#ifdef HAVE_UMFPACK
+          Matrix Control, Info;
+          void *Numeric = factorize (err, rcond, Control, Info,
+                                     sing_handler, calc_cond);
+
+          if (err == 0)
+            {
+              octave_idx_type b_nr = b.rows ();
+              octave_idx_type b_nc = b.cols ();
+              int status = 0;
+              double *control = Control.fortran_vec ();
+              double *info = Info.fortran_vec ();
+              const octave_idx_type *Ap = cidx ();
+              const octave_idx_type *Ai = ridx ();
+              const double *Ax = data ();
+
+              OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
+              OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
+
+              retval.resize (b_nr, b_nc);
+
+              OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
+              OCTAVE_LOCAL_BUFFER (double, Xz, b_nr);
+
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < b_nr; i++)
+                    {
+                      Complex c = b (i,j);
+                      Bx[i] = std::real (c);
+                      Bz[i] = std::imag (c);
+                    }
+
+                  status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap,
+                                             Ai, Ax, Xx, Bx, Numeric, control,
+                                             info);
+                  int status2 = UMFPACK_DNAME (solve) (UMFPACK_A,
+                                                  Ap, Ai, Ax, Xz, Bz, Numeric,
+                                                  control, info) ;
+
+                  if (status < 0 || status2 < 0)
+                    {
+                      (*current_liboctave_error_handler)
+                        ("SparseMatrix::solve solve failed");
+
+                      UMFPACK_DNAME (report_status) (control, status);
+
+                      err = -1;
+
+                      break;
+                    }
+
+                  for (octave_idx_type i = 0; i < b_nr; i++)
+                    retval(i, j) = Complex (Xx[i], Xz[i]);
+                }
+
+              UMFPACK_DNAME (report_info) (control, info);
+
+              UMFPACK_DNAME (free_numeric) (&Numeric);
+            }
+          else
+            mattype.mark_as_rectangular ();
+
+#else
+          (*current_liboctave_error_handler) ("UMFPACK not installed");
+#endif
+        }
+      else if (typ != MatrixType::Hermitian)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+SparseComplexMatrix
+SparseMatrix::fsolve (MatrixType &mattype, const SparseComplexMatrix& b,
+                      octave_idx_type& err, double& rcond,
+                      solve_singularity_handler sing_handler,
+                      bool calc_cond) const
+{
+  SparseComplexMatrix retval;
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  err = 0;
+
+  if (nr != nc || nr != b.rows ())
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch solution of linear equations");
+  else if (nr == 0 || b.cols () == 0)
+    retval = SparseComplexMatrix (nc, b.cols ());
+  else
+    {
+      // Print spparms("spumoni") info if requested
+      volatile int typ = mattype.type ();
+      mattype.info ();
+
+      if (typ == MatrixType::Hermitian)
+        {
+#ifdef HAVE_CHOLMOD
+          cholmod_common Common;
+          cholmod_common *cm = &Common;
+
+          // Setup initial parameters
+          CHOLMOD_NAME(start) (cm);
+          cm->prefer_zomplex = false;
+
+          double spu = octave_sparse_params::get_key ("spumoni");
+          if (spu == 0.)
+            {
+              cm->print = -1;
+              cm->print_function = 0;
+            }
+          else
+            {
+              cm->print = static_cast<int> (spu) + 2;
+              cm->print_function =&SparseCholPrint;
+            }
+
+          cm->error_handler = &SparseCholError;
+          cm->complex_divide = CHOLMOD_NAME(divcomplex);
+          cm->hypotenuse = CHOLMOD_NAME(hypot);
+
+          cm->final_ll = true;
+
+          cholmod_sparse Astore;
+          cholmod_sparse *A = &Astore;
+          double dummy;
+          A->nrow = nr;
+          A->ncol = nc;
+
+          A->p = cidx ();
+          A->i = ridx ();
+          A->nzmax = nnz ();
+          A->packed = true;
+          A->sorted = true;
+          A->nz = 0;
+#ifdef IDX_TYPE_LONG
+          A->itype = CHOLMOD_LONG;
+#else
+          A->itype = CHOLMOD_INT;
+#endif
+          A->dtype = CHOLMOD_DOUBLE;
+          A->stype = 1;
+          A->xtype = CHOLMOD_REAL;
+
+          if (nr < 1)
+            A->x = &dummy;
+          else
+            A->x = data ();
+
+          cholmod_sparse Bstore;
+          cholmod_sparse *B = &Bstore;
+          B->nrow = b.rows ();
+          B->ncol = b.cols ();
+          B->p = b.cidx ();
+          B->i = b.ridx ();
+          B->nzmax = b.nnz ();
+          B->packed = true;
+          B->sorted = true;
+          B->nz = 0;
+#ifdef IDX_TYPE_LONG
+          B->itype = CHOLMOD_LONG;
+#else
+          B->itype = CHOLMOD_INT;
+#endif
+          B->dtype = CHOLMOD_DOUBLE;
+          B->stype = 0;
+          B->xtype = CHOLMOD_COMPLEX;
+
+          if (b.rows () < 1 || b.cols () < 1)
+            B->x = &dummy;
+          else
+            B->x = b.data ();
+
+          cholmod_factor *L;
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          L = CHOLMOD_NAME(analyze) (A, cm);
+          CHOLMOD_NAME(factorize) (A, L, cm);
+          if (calc_cond)
+            rcond = CHOLMOD_NAME(rcond)(L, cm);
+          else
+            rcond = 1.0;
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+          if (rcond == 0.0)
+            {
+              // Either its indefinite or singular. Try UMFPACK
+              mattype.mark_as_unsymmetric ();
+              typ = MatrixType::Full;
+            }
+          else
+            {
+              volatile double rcond_plus_one = rcond + 1.0;
+
+              if (rcond_plus_one == 1.0 || xisnan (rcond))
+                {
+                  err = -2;
+
+                  if (sing_handler)
+                    {
+                      sing_handler (rcond);
+                      mattype.mark_as_rectangular ();
+                    }
+                  else
+                    (*current_liboctave_error_handler)
+                      ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
+                       rcond);
+
+                  return retval;
+                }
+
+              cholmod_sparse *X;
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+              retval = SparseComplexMatrix
+                (static_cast<octave_idx_type>(X->nrow),
+                 static_cast<octave_idx_type>(X->ncol),
+                 static_cast<octave_idx_type>(X->nzmax));
+              for (octave_idx_type j = 0;
+                   j <= static_cast<octave_idx_type>(X->ncol); j++)
+                retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
+              for (octave_idx_type j = 0;
+                   j < static_cast<octave_idx_type>(X->nzmax); j++)
+                {
+                  retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
+                  retval.xdata (j) = static_cast<Complex *>(X->x)[j];
+                }
+
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              CHOLMOD_NAME(free_sparse) (&X, cm);
+              CHOLMOD_NAME(free_factor) (&L, cm);
+              CHOLMOD_NAME(finish) (cm);
+              static char tmp[] = " ";
+              CHOLMOD_NAME(print_common) (tmp, cm);
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+            }
+#else
+          (*current_liboctave_warning_handler)
+            ("CHOLMOD not installed");
+
+          mattype.mark_as_unsymmetric ();
+          typ = MatrixType::Full;
+#endif
+        }
+
+      if (typ == MatrixType::Full)
+        {
+#ifdef HAVE_UMFPACK
+          Matrix Control, Info;
+          void *Numeric = factorize (err, rcond, Control, Info,
+                                     sing_handler, calc_cond);
+
+          if (err == 0)
+            {
+              octave_idx_type b_nr = b.rows ();
+              octave_idx_type b_nc = b.cols ();
+              int status = 0;
+              double *control = Control.fortran_vec ();
+              double *info = Info.fortran_vec ();
+              const octave_idx_type *Ap = cidx ();
+              const octave_idx_type *Ai = ridx ();
+              const double *Ax = data ();
+
+              OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
+              OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
+
+              // Take a first guess that the number of non-zero terms
+              // will be as many as in b
+              octave_idx_type x_nz = b.nnz ();
+              octave_idx_type ii = 0;
+              retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
+
+              OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
+              OCTAVE_LOCAL_BUFFER (double, Xz, b_nr);
+
+              retval.xcidx (0) = 0;
+              for (octave_idx_type j = 0; j < b_nc; j++)
+                {
+                  for (octave_idx_type i = 0; i < b_nr; i++)
+                    {
+                      Complex c = b (i,j);
+                      Bx[i] = std::real (c);
+                      Bz[i] = std::imag (c);
+                    }
+
+                  status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap,
+                                             Ai, Ax, Xx, Bx, Numeric, control,
+                                             info);
+                  int status2 = UMFPACK_DNAME (solve) (UMFPACK_A,
+                                                  Ap, Ai, Ax, Xz, Bz, Numeric,
+                                                  control, info) ;
+
+                  if (status < 0 || status2 < 0)
+                    {
+                      (*current_liboctave_error_handler)
+                        ("SparseMatrix::solve solve failed");
+
+                      UMFPACK_DNAME (report_status) (control, status);
+
+                      err = -1;
+
+                      break;
+                    }
+
+                  for (octave_idx_type i = 0; i < b_nr; i++)
+                    {
+                      Complex tmp = Complex (Xx[i], Xz[i]);
+                      if (tmp != 0.0)
+                        {
+                          if (ii == x_nz)
+                            {
+                              // Resize the sparse matrix
+                              octave_idx_type sz = x_nz * (b_nc - j) / b_nc;
+                              sz = (sz > 10 ? sz : 10) + x_nz;
+                              retval.change_capacity (sz);
+                              x_nz = sz;
+                            }
+                          retval.xdata (ii) = tmp;
+                          retval.xridx (ii++) = i;
+                        }
+                    }
+                  retval.xcidx (j+1) = ii;
+                }
+
+              retval.maybe_compress ();
+
+              UMFPACK_DNAME (report_info) (control, info);
+
+              UMFPACK_DNAME (free_numeric) (&Numeric);
+            }
+          else
+            mattype.mark_as_rectangular ();
+#else
+          (*current_liboctave_error_handler) ("UMFPACK not installed");
+#endif
+        }
+      else if (typ != MatrixType::Hermitian)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+    }
+
+  return retval;
+}
+
+Matrix
+SparseMatrix::solve (MatrixType &mattype, const Matrix& b) const
+{
+  octave_idx_type info;
+  double rcond;
+  return solve (mattype, b, info, rcond, 0);
+}
+
+Matrix
+SparseMatrix::solve (MatrixType &mattype, const Matrix& b,
+                     octave_idx_type& info) const
+{
+  double rcond;
+  return solve (mattype, b, info, rcond, 0);
+}
+
+Matrix
+SparseMatrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
+                     double& rcond) const
+{
+  return solve (mattype, b, info, rcond, 0);
+}
+
+Matrix
+SparseMatrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& err,
+                     double& rcond, solve_singularity_handler sing_handler,
+                     bool singular_fallback) const
+{
+  Matrix retval;
+  int typ = mattype.type (false);
+
+  if (typ == MatrixType::Unknown)
+    typ = mattype.type (*this);
+
+  // Only calculate the condition number for CHOLMOD/UMFPACK
+  if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
+    retval = dsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
+    retval = utsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
+    retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
+    retval = bsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Tridiagonal ||
+           typ == MatrixType::Tridiagonal_Hermitian)
+    retval = trisolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
+    retval = fsolve (mattype, b, err, rcond, sing_handler, true);
+  else if (typ != MatrixType::Rectangular)
+    {
+      (*current_liboctave_error_handler) ("unknown matrix type");
+      return Matrix ();
+    }
+
+  // Rectangular or one of the above solvers flags a singular matrix
+  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
+    {
+      rcond = 1.;
+#ifdef USE_QRSOLVE
+      retval = qrsolve (*this, b, err);
+#else
+      retval = dmsolve<Matrix, SparseMatrix, Matrix> (*this, b, err);
+#endif
+    }
+
+  return retval;
+}
+
+SparseMatrix
+SparseMatrix::solve (MatrixType &mattype, const SparseMatrix& b) const
+{
+  octave_idx_type info;
+  double rcond;
+  return solve (mattype, b, info, rcond, 0);
+}
+
+SparseMatrix
+SparseMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
+                     octave_idx_type& info) const
+{
+  double rcond;
+  return solve (mattype, b, info, rcond, 0);
+}
+
+SparseMatrix
+SparseMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
+                     octave_idx_type& info, double& rcond) const
+{
+  return solve (mattype, b, info, rcond, 0);
+}
+
+SparseMatrix
+SparseMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
+                     octave_idx_type& err, double& rcond,
+                     solve_singularity_handler sing_handler,
+                     bool singular_fallback) const
+{
+  SparseMatrix retval;
+  int typ = mattype.type (false);
+
+  if (typ == MatrixType::Unknown)
+    typ = mattype.type (*this);
+
+  if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
+    retval = dsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
+    retval = utsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
+    retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
+    retval = bsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Tridiagonal ||
+           typ == MatrixType::Tridiagonal_Hermitian)
+    retval = trisolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
+    retval = fsolve (mattype, b, err, rcond, sing_handler, true);
+  else if (typ != MatrixType::Rectangular)
+    {
+      (*current_liboctave_error_handler) ("unknown matrix type");
+      return SparseMatrix ();
+    }
+
+  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
+    {
+      rcond = 1.;
+#ifdef USE_QRSOLVE
+      retval = qrsolve (*this, b, err);
+#else
+      retval = dmsolve<SparseMatrix, SparseMatrix,
+        SparseMatrix> (*this, b, err);
+#endif
+    }
+
+  return retval;
+}
+
+ComplexMatrix
+SparseMatrix::solve (MatrixType &mattype, const ComplexMatrix& b) const
+{
+  octave_idx_type info;
+  double rcond;
+  return solve (mattype, b, info, rcond, 0);
+}
+
+ComplexMatrix
+SparseMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
+                            octave_idx_type& info) const
+{
+  double rcond;
+  return solve (mattype, b, info, rcond, 0);
+}
+
+ComplexMatrix
+SparseMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
+                     octave_idx_type& info, double& rcond) const
+{
+  return solve (mattype, b, info, rcond, 0);
+}
+
+ComplexMatrix
+SparseMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
+                     octave_idx_type& err, double& rcond,
+                     solve_singularity_handler sing_handler,
+                     bool singular_fallback) const
+{
+  ComplexMatrix retval;
+  int typ = mattype.type (false);
+
+  if (typ == MatrixType::Unknown)
+    typ = mattype.type (*this);
+
+  if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
+    retval = dsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
+    retval = utsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
+    retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
+    retval = bsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Tridiagonal ||
+           typ == MatrixType::Tridiagonal_Hermitian)
+    retval = trisolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
+    retval = fsolve (mattype, b, err, rcond, sing_handler, true);
+  else if (typ != MatrixType::Rectangular)
+    {
+      (*current_liboctave_error_handler) ("unknown matrix type");
+      return ComplexMatrix ();
+    }
+
+  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
+    {
+      rcond = 1.;
+#ifdef USE_QRSOLVE
+      retval = qrsolve (*this, b, err);
+#else
+      retval = dmsolve<ComplexMatrix, SparseMatrix,
+        ComplexMatrix> (*this, b, err);
+#endif
+    }
+
+  return retval;
+}
+
+SparseComplexMatrix
+SparseMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b) const
+{
+  octave_idx_type info;
+  double rcond;
+  return solve (mattype, b, info, rcond, 0);
+}
+
+SparseComplexMatrix
+SparseMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b,
+                     octave_idx_type& info) const
+{
+  double rcond;
+  return solve (mattype, b, info, rcond, 0);
+}
+
+SparseComplexMatrix
+SparseMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b,
+                     octave_idx_type& info, double& rcond) const
+{
+  return solve (mattype, b, info, rcond, 0);
+}
+
+SparseComplexMatrix
+SparseMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b,
+                     octave_idx_type& err, double& rcond,
+                     solve_singularity_handler sing_handler,
+                     bool singular_fallback) const
+{
+  SparseComplexMatrix retval;
+  int typ = mattype.type (false);
+
+  if (typ == MatrixType::Unknown)
+    typ = mattype.type (*this);
+
+  if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
+    retval = dsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
+    retval = utsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
+    retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
+    retval = bsolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Tridiagonal ||
+           typ == MatrixType::Tridiagonal_Hermitian)
+    retval = trisolve (mattype, b, err, rcond, sing_handler, false);
+  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
+    retval = fsolve (mattype, b, err, rcond, sing_handler, true);
+  else if (typ != MatrixType::Rectangular)
+    {
+      (*current_liboctave_error_handler) ("unknown matrix type");
+      return SparseComplexMatrix ();
+    }
+
+  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
+    {
+      rcond = 1.;
+#ifdef USE_QRSOLVE
+      retval = qrsolve (*this, b, err);
+#else
+      retval = dmsolve<SparseComplexMatrix, SparseMatrix,
+        SparseComplexMatrix> (*this, b, err);
+#endif
+    }
+
+  return retval;
+}
+
+ColumnVector
+SparseMatrix::solve (MatrixType &mattype, const ColumnVector& b) const
+{
+  octave_idx_type info; double rcond;
+  return solve (mattype, b, info, rcond);
+}
+
+ColumnVector
+SparseMatrix::solve (MatrixType &mattype, const ColumnVector& b, octave_idx_type& info) const
+{
+  double rcond;
+  return solve (mattype, b, info, rcond);
+}
+
+ColumnVector
+SparseMatrix::solve (MatrixType &mattype, const ColumnVector& b, octave_idx_type& info, double& rcond) const
+{
+  return solve (mattype, b, info, rcond, 0);
+}
+
+ColumnVector
+SparseMatrix::solve (MatrixType &mattype, const ColumnVector& b, octave_idx_type& info, double& rcond,
+               solve_singularity_handler sing_handler) const
+{
+  Matrix tmp (b);
+  return solve (mattype, tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0));
+}
+
+ComplexColumnVector
+SparseMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b) const
+{
+  octave_idx_type info;
+  double rcond;
+  return solve (mattype, b, info, rcond, 0);
+}
+
+ComplexColumnVector
+SparseMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b, octave_idx_type& info) const
+{
+  double rcond;
+  return solve (mattype, b, info, rcond, 0);
+}
+
+ComplexColumnVector
+SparseMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b, octave_idx_type& info,
+                     double& rcond) const
+{
+  return solve (mattype, b, info, rcond, 0);
+}
+
+ComplexColumnVector
+SparseMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b, octave_idx_type& info, double& rcond,
+               solve_singularity_handler sing_handler) const
+{
+  ComplexMatrix tmp (b);
+  return solve (mattype, tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0));
+}
+
+Matrix
+SparseMatrix::solve (const Matrix& b) const
+{
+  octave_idx_type info;
+  double rcond;
+  return solve (b, info, rcond, 0);
+}
+
+Matrix
+SparseMatrix::solve (const Matrix& b, octave_idx_type& info) const
+{
+  double rcond;
+  return solve (b, info, rcond, 0);
+}
+
+Matrix
+SparseMatrix::solve (const Matrix& b, octave_idx_type& info,
+                     double& rcond) const
+{
+  return solve (b, info, rcond, 0);
+}
+
+Matrix
+SparseMatrix::solve (const Matrix& b, octave_idx_type& err,
+                     double& rcond,
+                     solve_singularity_handler sing_handler) const
+{
+  MatrixType mattype (*this);
+  return solve (mattype, b, err, rcond, sing_handler);
+}
+
+SparseMatrix
+SparseMatrix::solve (const SparseMatrix& b) const
+{
+  octave_idx_type info;
+  double rcond;
+  return solve (b, info, rcond, 0);
+}
+
+SparseMatrix
+SparseMatrix::solve (const SparseMatrix& b,
+                     octave_idx_type& info) const
+{
+  double rcond;
+  return solve (b, info, rcond, 0);
+}
+
+SparseMatrix
+SparseMatrix::solve (const SparseMatrix& b,
+                     octave_idx_type& info, double& rcond) const
+{
+  return solve (b, info, rcond, 0);
+}
+
+SparseMatrix
+SparseMatrix::solve (const SparseMatrix& b,
+                     octave_idx_type& err, double& rcond,
+                     solve_singularity_handler sing_handler) const
+{
+  MatrixType mattype (*this);
+  return solve (mattype, b, err, rcond, sing_handler);
+}
+
+ComplexMatrix
+SparseMatrix::solve (const ComplexMatrix& b,
+                            octave_idx_type& info) const
+{
+  double rcond;
+  return solve (b, info, rcond, 0);
+}
+
+ComplexMatrix
+SparseMatrix::solve (const ComplexMatrix& b,
+                     octave_idx_type& info, double& rcond) const
+{
+  return solve (b, info, rcond, 0);
+}
+
+ComplexMatrix
+SparseMatrix::solve (const ComplexMatrix& b,
+                     octave_idx_type& err, double& rcond,
+                     solve_singularity_handler sing_handler) const
+{
+  MatrixType mattype (*this);
+  return solve (mattype, b, err, rcond, sing_handler);
+}
+
+SparseComplexMatrix
+SparseMatrix::solve (const SparseComplexMatrix& b) const
+{
+  octave_idx_type info;
+  double rcond;
+  return solve (b, info, rcond, 0);
+}
+
+SparseComplexMatrix
+SparseMatrix::solve (const SparseComplexMatrix& b,
+                     octave_idx_type& info) const
+{
+  double rcond;
+  return solve (b, info, rcond, 0);
+}
+
+SparseComplexMatrix
+SparseMatrix::solve (const SparseComplexMatrix& b,
+                     octave_idx_type& info, double& rcond) const
+{
+  return solve (b, info, rcond, 0);
+}
+
+SparseComplexMatrix
+SparseMatrix::solve (const SparseComplexMatrix& b,
+                     octave_idx_type& err, double& rcond,
+                     solve_singularity_handler sing_handler) const
+{
+  MatrixType mattype (*this);
+  return solve (mattype, b, err, rcond, sing_handler);
+}
+
+ColumnVector
+SparseMatrix::solve (const ColumnVector& b) const
+{
+  octave_idx_type info; double rcond;
+  return solve (b, info, rcond);
+}
+
+ColumnVector
+SparseMatrix::solve (const ColumnVector& b, octave_idx_type& info) const
+{
+  double rcond;
+  return solve (b, info, rcond);
+}
+
+ColumnVector
+SparseMatrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond) const
+{
+  return solve (b, info, rcond, 0);
+}
+
+ColumnVector
+SparseMatrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond,
+               solve_singularity_handler sing_handler) const
+{
+  Matrix tmp (b);
+  return solve (tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0));
+}
+
+ComplexColumnVector
+SparseMatrix::solve (const ComplexColumnVector& b) const
+{
+  octave_idx_type info;
+  double rcond;
+  return solve (b, info, rcond, 0);
+}
+
+ComplexColumnVector
+SparseMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info) const
+{
+  double rcond;
+  return solve (b, info, rcond, 0);
+}
+
+ComplexColumnVector
+SparseMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
+                     double& rcond) const
+{
+  return solve (b, info, rcond, 0);
+}
+
+ComplexColumnVector
+SparseMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcond,
+               solve_singularity_handler sing_handler) const
+{
+  ComplexMatrix tmp (b);
+  return solve (tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0));
+}
+
+// other operations.
+
+bool
+SparseMatrix::any_element_is_negative (bool neg_zero) const
+{
+  octave_idx_type nel = nnz ();
+
+  if (neg_zero)
+    {
+      for (octave_idx_type i = 0; i < nel; i++)
+        if (lo_ieee_signbit (data (i)))
+          return true;
+    }
+  else
+    {
+      for (octave_idx_type i = 0; i < nel; i++)
+        if (data (i) < 0)
+          return true;
+    }
+
+  return false;
+}
+
+bool
+SparseMatrix::any_element_is_nan (void) const
+{
+  octave_idx_type nel = nnz ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    {
+      double val = data (i);
+      if (xisnan (val))
+        return true;
+    }
+
+  return false;
+}
+
+bool
+SparseMatrix::any_element_is_inf_or_nan (void) const
+{
+  octave_idx_type nel = nnz ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    {
+      double val = data (i);
+      if (xisinf (val) || xisnan (val))
+        return true;
+    }
+
+  return false;
+}
+
+bool
+SparseMatrix::any_element_not_one_or_zero (void) const
+{
+  octave_idx_type nel = nnz ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    {
+      double val = data (i);
+      if (val != 0.0 && val != 1.0)
+        return true;
+    }
+
+  return false;
+}
+
+bool
+SparseMatrix::all_elements_are_zero (void) const
+{
+  octave_idx_type nel = nnz ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    if (data (i) != 0)
+      return false;
+
+  return true;
+}
+
+bool
+SparseMatrix::all_elements_are_int_or_inf_or_nan (void) const
+{
+  octave_idx_type nel = nnz ();
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    {
+      double val = data (i);
+      if (xisnan (val) || D_NINT (val) == val)
+        continue;
+      else
+        return false;
+    }
+
+  return true;
+}
+
+// Return nonzero if any element of M is not an integer.  Also extract
+// the largest and smallest values and return them in MAX_VAL and MIN_VAL.
+
+bool
+SparseMatrix::all_integers (double& max_val, double& min_val) const
+{
+  octave_idx_type nel = nnz ();
+
+  if (nel == 0)
+    return false;
+
+  max_val = data (0);
+  min_val = data (0);
+
+  for (octave_idx_type i = 0; i < nel; i++)
+    {
+      double val = data (i);
+
+      if (val > max_val)
+        max_val = val;
+
+      if (val < min_val)
+        min_val = val;
+
+      if (D_NINT (val) != val)
+        return false;
+    }
+
+  return true;
+}
+
+bool
+SparseMatrix::too_large_for_float (void) const
+{
+  return test_any (xtoo_large_for_float);
+}
+
+SparseBoolMatrix
+SparseMatrix::operator ! (void) const
+{
+  if (any_element_is_nan ())
+    gripe_nan_to_logical_conversion ();
+
+  octave_idx_type nr = rows ();
+  octave_idx_type nc = cols ();
+  octave_idx_type nz1 = nnz ();
+  octave_idx_type nz2 = nr*nc - nz1;
+
+  SparseBoolMatrix r (nr, nc, nz2);
+
+  octave_idx_type ii = 0;
+  octave_idx_type jj = 0;
+  r.cidx (0) = 0;
+  for (octave_idx_type i = 0; i < nc; i++)
+    {
+      for (octave_idx_type j = 0; j < nr; j++)
+        {
+          if (jj < cidx (i+1) && ridx (jj) == j)
+            jj++;
+          else
+            {
+              r.data (ii) = true;
+              r.ridx (ii++) = j;
+            }
+        }
+      r.cidx (i+1) = ii;
+    }
+
+  return r;
+}
+
+// FIXME Do these really belong here?  Maybe they should be
+// in a base class?
+
+SparseBoolMatrix
+SparseMatrix::all (int dim) const
+{
+  SPARSE_ALL_OP (dim);
+}
+
+SparseBoolMatrix
+SparseMatrix::any (int dim) const
+{
+  SPARSE_ANY_OP (dim);
+}
+
+SparseMatrix
+SparseMatrix::cumprod (int dim) const
+{
+  SPARSE_CUMPROD (SparseMatrix, double, cumprod);
+}
+
+SparseMatrix
+SparseMatrix::cumsum (int dim) const
+{
+  SPARSE_CUMSUM (SparseMatrix, double, cumsum);
+}
+
+SparseMatrix
+SparseMatrix::prod (int dim) const
+{
+  if ((rows () == 1 && dim == -1) || dim == 1)
+    return transpose (). prod (0). transpose ();
+  else
+    {
+      SPARSE_REDUCTION_OP (SparseMatrix, double, *=,
+                           (cidx (j+1) - cidx (j) < nr ? 0.0 : 1.0), 1.0);
+    }
+}
+
+SparseMatrix
+SparseMatrix::sum (int dim) const
+{
+  SPARSE_REDUCTION_OP (SparseMatrix, double, +=, 0.0, 0.0);
+}
+
+SparseMatrix
+SparseMatrix::sumsq (int dim) const
+{
+#define ROW_EXPR \
+  double d = data (i); \
+  tmp[ridx (i)] += d * d
+
+#define COL_EXPR \
+  double d = data (i); \
+  tmp[j] += d * d
+
+  SPARSE_BASE_REDUCTION_OP (SparseMatrix, double, ROW_EXPR, COL_EXPR,
+                            0.0, 0.0);
+
+#undef ROW_EXPR
+#undef COL_EXPR
+}
+
+SparseMatrix
+SparseMatrix::abs (void) const
+{
+  octave_idx_type nz = nnz ();
+
+  SparseMatrix retval (*this);
+
+  for (octave_idx_type i = 0; i < nz; i++)
+    retval.data (i) = fabs (retval.data (i));
+
+  return retval;
+}
+
+SparseMatrix
+SparseMatrix::diag (octave_idx_type k) const
+{
+  return MSparse<double>::diag (k);
+}
+
+Matrix
+SparseMatrix::matrix_value (void) const
+{
+  return Sparse<double>::array_value ();
+}
+
+std::ostream&
+operator << (std::ostream& os, const SparseMatrix& a)
+{
+  octave_idx_type nc = a.cols ();
+
+   // add one to the printed indices to go from
+   //  zero-based to one-based arrays
+   for (octave_idx_type j = 0; j < nc; j++)
+     {
+       octave_quit ();
+       for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
+         {
+           os << a.ridx (i) + 1 << " "  << j + 1 << " ";
+           octave_write_double (os, a.data (i));
+           os << "\n";
+         }
+     }
+
+  return os;
+}
+
+std::istream&
+operator >> (std::istream& is, SparseMatrix& a)
+{
+  typedef SparseMatrix::element_type elt_type;
+
+  return read_sparse_matrix<elt_type> (is, a, octave_read_value<double>);
+}
+
+SparseMatrix
+SparseMatrix::squeeze (void) const
+{
+  return MSparse<double>::squeeze ();
+}
+
+SparseMatrix
+SparseMatrix::reshape (const dim_vector& new_dims) const
+{
+  return MSparse<double>::reshape (new_dims);
+}
+
+SparseMatrix
+SparseMatrix::permute (const Array<octave_idx_type>& vec, bool inv) const
+{
+  return MSparse<double>::permute (vec, inv);
+}
+
+SparseMatrix
+SparseMatrix::ipermute (const Array<octave_idx_type>& vec) const
+{
+  return MSparse<double>::ipermute (vec);
+}
+
+// matrix by matrix -> matrix operations
+
+SparseMatrix
+operator * (const SparseMatrix& m, const SparseMatrix& a)
+{
+  SPARSE_SPARSE_MUL (SparseMatrix, double, double);
+}
+
+Matrix
+operator * (const Matrix& m, const SparseMatrix& a)
+{
+  FULL_SPARSE_MUL (Matrix, double, 0.);
+}
+
+Matrix
+mul_trans (const Matrix& m, const SparseMatrix& a)
+{
+  FULL_SPARSE_MUL_TRANS (Matrix, double, 0., );
+}
+
+Matrix
+operator * (const SparseMatrix& m, const Matrix& a)
+{
+  SPARSE_FULL_MUL (Matrix, double, 0.);
+}
+
+Matrix
+trans_mul (const SparseMatrix& m, const Matrix& a)
+{
+  SPARSE_FULL_TRANS_MUL (Matrix, double, 0., );
+}
+
+// diag * sparse and sparse * diag
+
+SparseMatrix
+operator * (const DiagMatrix& d, const SparseMatrix& a)
+{
+  return do_mul_dm_sm<SparseMatrix> (d, a);
+}
+
+SparseMatrix
+operator * (const SparseMatrix& a, const DiagMatrix& d)
+{
+  return do_mul_sm_dm<SparseMatrix> (a, d);
+}
+
+SparseMatrix
+operator + (const DiagMatrix& d, const SparseMatrix& a)
+{
+  return do_add_dm_sm<SparseMatrix> (d, a);
+}
+
+SparseMatrix
+operator - (const DiagMatrix& d, const SparseMatrix& a)
+{
+  return do_sub_dm_sm<SparseMatrix> (d, a);
+}
+
+SparseMatrix
+operator + (const SparseMatrix& a, const DiagMatrix& d)
+{
+  return do_add_sm_dm<SparseMatrix> (a, d);
+}
+
+SparseMatrix
+operator - (const SparseMatrix& a, const DiagMatrix& d)
+{
+  return do_sub_sm_dm<SparseMatrix> (a, d);
+}
+
+// perm * sparse and sparse * perm
+
+SparseMatrix
+operator * (const PermMatrix& p, const SparseMatrix& a)
+{
+  return octinternal_do_mul_pm_sm (p, a);
+}
+
+SparseMatrix
+operator * (const SparseMatrix& a, const PermMatrix& p)
+{
+  return octinternal_do_mul_sm_pm (a, p);
+}
+
+// FIXME -- it would be nice to share code among the min/max
+// functions below.
+
+#define EMPTY_RETURN_CHECK(T) \
+  if (nr == 0 || nc == 0) \
+    return T (nr, nc);
+
+SparseMatrix
+min (double d, const SparseMatrix& m)
+{
+  SparseMatrix result;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (SparseMatrix);
+
+  // Count the number of non-zero elements
+  if (d < 0.)
+    {
+      result = SparseMatrix (nr, nc, d);
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
+          {
+            double tmp = xmin (d, m.data (i));
+            if (tmp != 0.)
+              {
+                octave_idx_type idx = m.ridx (i) + j * nr;
+                result.xdata (idx) = tmp;
+                result.xridx (idx) = m.ridx (i);
+              }
+          }
+    }
+  else
+    {
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
+          if (xmin (d, m.data (i)) != 0.)
+            nel++;
+
+      result = SparseMatrix (nr, nc, nel);
+
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
+            {
+              double tmp = xmin (d, m.data (i));
+
+              if (tmp != 0.)
+                {
+                  result.xdata (ii) = tmp;
+                  result.xridx (ii++) = m.ridx (i);
+                }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+
+  return result;
+}
+
+SparseMatrix
+min (const SparseMatrix& m, double d)
+{
+  return min (d, m);
+}
+
+SparseMatrix
+min (const SparseMatrix& a, const SparseMatrix& b)
+{
+  SparseMatrix r;
+
+  if ((a.rows () == b.rows ()) && (a.cols () == b.cols ()))
+    {
+      octave_idx_type a_nr = a.rows ();
+      octave_idx_type a_nc = a.cols ();
+
+      octave_idx_type b_nr = b.rows ();
+      octave_idx_type b_nc = b.cols ();
+
+      if (a_nr != b_nr || a_nc != b_nc)
+        gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
+      else
+        {
+          r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
+
+          octave_idx_type jx = 0;
+          r.cidx (0) = 0;
+          for (octave_idx_type i = 0 ; i < a_nc ; i++)
+            {
+              octave_idx_type  ja = a.cidx (i);
+              octave_idx_type  ja_max = a.cidx (i+1);
+              bool ja_lt_max= ja < ja_max;
+
+              octave_idx_type  jb = b.cidx (i);
+              octave_idx_type  jb_max = b.cidx (i+1);
+              bool jb_lt_max = jb < jb_max;
+
+              while (ja_lt_max || jb_lt_max )
+                {
+                  octave_quit ();
+                  if ((! jb_lt_max) ||
+                      (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
+                    {
+                      double tmp = xmin (a.data (ja), 0.);
+                      if (tmp != 0.)
+                        {
+                          r.ridx (jx) = a.ridx (ja);
+                          r.data (jx) = tmp;
+                          jx++;
+                        }
+                      ja++;
+                      ja_lt_max= ja < ja_max;
+                    }
+                  else if (( !ja_lt_max ) ||
+                           (jb_lt_max && (b.ridx (jb) < a.ridx (ja)) ) )
+                    {
+                      double tmp = xmin (0., b.data (jb));
+                      if (tmp != 0.)
+                        {
+                          r.ridx (jx) = b.ridx (jb);
+                          r.data (jx) = tmp;
+                          jx++;
+                        }
+                      jb++;
+                      jb_lt_max= jb < jb_max;
+                    }
+                  else
+                    {
+                      double tmp = xmin (a.data (ja), b.data (jb));
+                      if (tmp != 0.)
+                        {
+                          r.data (jx) = tmp;
+                          r.ridx (jx) = a.ridx (ja);
+                          jx++;
+                        }
+                      ja++;
+                      ja_lt_max= ja < ja_max;
+                      jb++;
+                      jb_lt_max= jb < jb_max;
+                    }
+                }
+              r.cidx (i+1) = jx;
+            }
+
+          r.maybe_compress ();
+        }
+    }
+  else
+    (*current_liboctave_error_handler) ("matrix size mismatch");
+
+  return r;
+}
+
+SparseMatrix
+max (double d, const SparseMatrix& m)
+{
+  SparseMatrix result;
+
+  octave_idx_type nr = m.rows ();
+  octave_idx_type nc = m.columns ();
+
+  EMPTY_RETURN_CHECK (SparseMatrix);
+
+  // Count the number of non-zero elements
+  if (d > 0.)
+    {
+      result = SparseMatrix (nr, nc, d);
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
+          {
+            double tmp = xmax (d, m.data (i));
+
+            if (tmp != 0.)
+              {
+                octave_idx_type idx = m.ridx (i) + j * nr;
+                result.xdata (idx) = tmp;
+                result.xridx (idx) = m.ridx (i);
+              }
+          }
+    }
+  else
+    {
+      octave_idx_type nel = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
+          if (xmax (d, m.data (i)) != 0.)
+            nel++;
+
+      result = SparseMatrix (nr, nc, nel);
+
+      octave_idx_type ii = 0;
+      result.xcidx (0) = 0;
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
+            {
+              double tmp = xmax (d, m.data (i));
+              if (tmp != 0.)
+                {
+                  result.xdata (ii) = tmp;
+                  result.xridx (ii++) = m.ridx (i);
+                }
+            }
+          result.xcidx (j+1) = ii;
+        }
+    }
+
+  return result;
+}
+
+SparseMatrix
+max (const SparseMatrix& m, double d)
+{
+  return max (d, m);
+}
+
+SparseMatrix
+max (const SparseMatrix& a, const SparseMatrix& b)
+{
+  SparseMatrix r;
+
+  if ((a.rows () == b.rows ()) && (a.cols () == b.cols ()))
+    {
+      octave_idx_type a_nr = a.rows ();
+      octave_idx_type a_nc = a.cols ();
+
+      octave_idx_type b_nr = b.rows ();
+      octave_idx_type b_nc = b.cols ();
+
+      if (a_nr != b_nr || a_nc != b_nc)
+        gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
+      else
+        {
+          r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
+
+          octave_idx_type jx = 0;
+          r.cidx (0) = 0;
+          for (octave_idx_type i = 0 ; i < a_nc ; i++)
+            {
+              octave_idx_type  ja = a.cidx (i);
+              octave_idx_type  ja_max = a.cidx (i+1);
+              bool ja_lt_max= ja < ja_max;
+
+              octave_idx_type  jb = b.cidx (i);
+              octave_idx_type  jb_max = b.cidx (i+1);
+              bool jb_lt_max = jb < jb_max;
+
+              while (ja_lt_max || jb_lt_max )
+                {
+                  octave_quit ();
+                  if ((! jb_lt_max) ||
+                      (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
+                    {
+                      double tmp = xmax (a.data (ja), 0.);
+                      if (tmp != 0.)
+                        {
+                          r.ridx (jx) = a.ridx (ja);
+                          r.data (jx) = tmp;
+                          jx++;
+                        }
+                      ja++;
+                      ja_lt_max= ja < ja_max;
+                    }
+                  else if (( !ja_lt_max ) ||
+                           (jb_lt_max && (b.ridx (jb) < a.ridx (ja)) ) )
+                    {
+                      double tmp = xmax (0., b.data (jb));
+                      if (tmp != 0.)
+                        {
+                          r.ridx (jx) = b.ridx (jb);
+                          r.data (jx) = tmp;
+                          jx++;
+                        }
+                      jb++;
+                      jb_lt_max= jb < jb_max;
+                    }
+                  else
+                    {
+                      double tmp = xmax (a.data (ja), b.data (jb));
+                      if (tmp != 0.)
+                        {
+                          r.data (jx) = tmp;
+                          r.ridx (jx) = a.ridx (ja);
+                          jx++;
+                        }
+                      ja++;
+                      ja_lt_max= ja < ja_max;
+                      jb++;
+                      jb_lt_max= jb < jb_max;
+                    }
+                }
+              r.cidx (i+1) = jx;
+            }
+
+          r.maybe_compress ();
+        }
+    }
+  else
+    (*current_liboctave_error_handler) ("matrix size mismatch");
+
+  return r;
+}
+
+SPARSE_SMS_CMP_OPS (SparseMatrix, 0.0, , double, 0.0, )
+SPARSE_SMS_BOOL_OPS (SparseMatrix, double, 0.0)
+
+SPARSE_SSM_CMP_OPS (double, 0.0, , SparseMatrix, 0.0, )
+SPARSE_SSM_BOOL_OPS (double, SparseMatrix, 0.0)
+
+SPARSE_SMSM_CMP_OPS (SparseMatrix, 0.0, , SparseMatrix, 0.0, )
+SPARSE_SMSM_BOOL_OPS (SparseMatrix, SparseMatrix, 0.0)