--- a/liboctave/CSparse.cc
+++ b/liboctave/CSparse.cc
@@ -47,6 +47,13 @@
 
 #include "oct-sort.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"
 {
@@ -82,7 +89,7 @@
   F77_RET_T
   F77_FUNC (zpbcon, ZPBCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&, 
 			     const octave_idx_type&, Complex*, const octave_idx_type&, 
-			     const double&, double&, Complex*, octave_idx_type*, octave_idx_type&
+			     const double&, double&, Complex*, double*, octave_idx_type&
 			     F77_CHAR_ARG_LEN_DECL);
 
   F77_RET_T
@@ -106,33 +113,33 @@
 }
 
 SparseComplexMatrix::SparseComplexMatrix (const SparseMatrix& a)
-  : MSparse<Complex> (a.rows (), a.cols (), a.nzmax ())
+  : MSparse<Complex> (a.rows (), a.cols (), a.nnz ())
 {
   octave_idx_type nc = cols ();
-  octave_idx_type nz = nzmax ();
+  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);
+      data (i) = Complex (a.data (i));
       ridx (i) = a.ridx (i);
     }
 }
 
 SparseComplexMatrix::SparseComplexMatrix (const SparseBoolMatrix& a)
-  : MSparse<Complex> (a.rows (), a.cols (), a.nzmax ())
+  : MSparse<Complex> (a.rows (), a.cols (), a.nnz ())
 {
   octave_idx_type nc = cols ();
-  octave_idx_type nz = nzmax ();
+  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);
+      data (i) = Complex (a.data (i));
       ridx (i) = a.ridx (i);
     }
 }
@@ -142,10 +149,10 @@
 {
   octave_idx_type nr = rows ();
   octave_idx_type nc = cols ();
-  octave_idx_type nz = nzmax ();
+  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.nzmax ();
+  octave_idx_type nz_a = a.nnz ();
 
   if (nr != nr_a || nc != nc_a || nz != nz_a)
     return false;
@@ -546,7 +553,7 @@
 {
   octave_idx_type nr = rows ();
   octave_idx_type nc = cols ();
-  octave_idx_type nz = nzmax ();
+  octave_idx_type nz = nnz ();
   SparseComplexMatrix retval (nc, nr, nz);
 
   OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr + 1);
@@ -580,7 +587,7 @@
 {
   octave_idx_type nr = a.rows ();
   octave_idx_type nc = a.cols ();
-  octave_idx_type nz = a.nzmax ();
+  octave_idx_type nz = a.nnz ();
   SparseComplexMatrix retval (nc, nr, nz);
 
   for (octave_idx_type i = 0; i < nc + 1; i++)
@@ -713,7 +720,7 @@
 
 	  if (typ == SparseType::Upper || typ == SparseType::Lower)
 	    {
-	      octave_idx_type nz = nzmax ();
+	      octave_idx_type nz = nnz ();
 	      octave_idx_type cx = 0;
 	      octave_idx_type nz2 = nz;
 	      retval = SparseComplexMatrix (nr, nc, nz2);
@@ -798,7 +805,7 @@
 	    }
 	  else
 	    {
-	      octave_idx_type nz = nzmax ();
+	      octave_idx_type nz = nnz ();
 	      octave_idx_type cx = 0;
 	      octave_idx_type nz2 = nz;
 	      retval = SparseComplexMatrix (nr, nc, nz2);
@@ -1118,8 +1125,8 @@
 
 ComplexMatrix
 SparseComplexMatrix::dsolve (SparseType &mattype, const Matrix& b,
-			     octave_idx_type& err, 
-			     double& rcond, solve_singularity_handler) const
+			     octave_idx_type& err, double& rcond, 
+			     solve_singularity_handler, bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -1151,16 +1158,21 @@
 		for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
 		  retval(k,j) = b(ridx(i),j) / data (i);
 	    
-	  double dmax = 0., dmin = octave_Inf; 
-	  for (octave_idx_type i = 0; i < nm; i++)
-	    {
-	      double tmp = std::abs(data(i));
-	      if (tmp > dmax)
-		dmax = tmp;
-	      if (tmp < dmin)
-		dmin = tmp;
-	    }
-	  rcond = dmin / dmax;
+	  if (calc_cond)
+	    {
+	      double dmax = 0., dmin = octave_Inf; 
+	      for (octave_idx_type i = 0; i < nm; i++)
+		{
+		  double tmp = std::abs(data(i));
+		  if (tmp > dmax)
+		    dmax = tmp;
+		  if (tmp < dmin)
+		    dmin = tmp;
+		}
+	      rcond = dmin / dmax;
+	    }
+	  else
+	    rcond = 1.0;
 	}
       else
 	(*current_liboctave_error_handler) ("incorrect matrix type");
@@ -1172,7 +1184,8 @@
 SparseComplexMatrix
 SparseComplexMatrix::dsolve (SparseType &mattype, const SparseMatrix& b,
 			     octave_idx_type& err, double& rcond, 
-			     solve_singularity_handler) const
+			     solve_singularity_handler,
+			     bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -1194,7 +1207,7 @@
 	  typ == SparseType::Permuted_Diagonal)
 	{
 	  octave_idx_type b_nc = b.cols ();
-	  octave_idx_type b_nz = b.nzmax ();
+	  octave_idx_type b_nz = b.nnz ();
 	  retval = SparseComplexMatrix (nc, b_nc, b_nz);
 
 	  retval.xcidx(0) = 0;
@@ -1204,6 +1217,8 @@
 	      {
 		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));
 		  }
@@ -1232,16 +1247,21 @@
 		retval.xcidx(j+1) = ii;
 	      }
 	    
-	  double dmax = 0., dmin = octave_Inf; 
-	  for (octave_idx_type i = 0; i < nm; i++)
-	    {
-	      double tmp = std::abs(data(i));
-	      if (tmp > dmax)
-		dmax = tmp;
-	      if (tmp < dmin)
-		dmin = tmp;
-	    }
-	  rcond = dmin / dmax;
+	  if (calc_cond)
+	    {
+	      double dmax = 0., dmin = octave_Inf; 
+	      for (octave_idx_type i = 0; i < nm; i++)
+		{
+		  double tmp = std::abs(data(i));
+		  if (tmp > dmax)
+		    dmax = tmp;
+		  if (tmp < dmin)
+		    dmin = tmp;
+		}
+	      rcond = dmin / dmax;
+	    }
+	  else
+	    rcond = 1.0;
 	}
       else
 	(*current_liboctave_error_handler) ("incorrect matrix type");
@@ -1253,7 +1273,8 @@
 ComplexMatrix
 SparseComplexMatrix::dsolve (SparseType &mattype, const ComplexMatrix& b,
 			     octave_idx_type& err, double& rcond, 
-			     solve_singularity_handler) const
+			     solve_singularity_handler,
+			     bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -1285,16 +1306,21 @@
 		for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
 		  retval(k,j) = b(ridx(i),j) / data (i);
 	    
-	  double dmax = 0., dmin = octave_Inf; 
-	  for (octave_idx_type i = 0; i < nr; i++)
-	    {
-	      double tmp = std::abs(data(i));
-	      if (tmp > dmax)
-		dmax = tmp;
-	      if (tmp < dmin)
-		dmin = tmp;
-	    }
-	  rcond = dmin / dmax;
+	  if (calc_cond)
+	    {
+	      double dmax = 0., dmin = octave_Inf; 
+	      for (octave_idx_type i = 0; i < nr; i++)
+		{
+		  double tmp = std::abs(data(i));
+		  if (tmp > dmax)
+		    dmax = tmp;
+		  if (tmp < dmin)
+		    dmin = tmp;
+		}
+	      rcond = dmin / dmax;
+	    }
+	  else
+	    rcond = 1.0;
 	}
       else
 	(*current_liboctave_error_handler) ("incorrect matrix type");
@@ -1306,7 +1332,8 @@
 SparseComplexMatrix
 SparseComplexMatrix::dsolve (SparseType &mattype, const SparseComplexMatrix& b,
 			     octave_idx_type& err, double& rcond, 
-			     solve_singularity_handler) const
+			     solve_singularity_handler,
+			     bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -1328,7 +1355,7 @@
 	  typ == SparseType::Permuted_Diagonal)
 	{
 	  octave_idx_type b_nc = b.cols ();
-	  octave_idx_type b_nz = b.nzmax ();
+	  octave_idx_type b_nz = b.nnz ();
 	  retval = SparseComplexMatrix (nc, b_nc, b_nz);
 
 	  retval.xcidx(0) = 0;
@@ -1338,6 +1365,8 @@
 	      {
 		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));
 		  }
@@ -1366,16 +1395,21 @@
 		retval.xcidx(j+1) = ii;
 	      }
 	    
-	  double dmax = 0., dmin = octave_Inf; 
-	  for (octave_idx_type i = 0; i < nm; i++)
-	    {
-	      double tmp = std::abs(data(i));
-	      if (tmp > dmax)
-		dmax = tmp;
-	      if (tmp < dmin)
-		dmin = tmp;
-	    }
-	  rcond = dmin / dmax;
+	  if (calc_cond)
+	    {
+	      double dmax = 0., dmin = octave_Inf; 
+	      for (octave_idx_type i = 0; i < nm; i++)
+		{
+		  double tmp = std::abs(data(i));
+		  if (tmp > dmax)
+		    dmax = tmp;
+		  if (tmp < dmin)
+		    dmin = tmp;
+		}
+	      rcond = dmin / dmax;
+	    }
+	  else
+	    rcond = 1.0;
 	}
       else
 	(*current_liboctave_error_handler) ("incorrect matrix type");
@@ -1387,7 +1421,8 @@
 ComplexMatrix
 SparseComplexMatrix::utsolve (SparseType &mattype, const Matrix& b,
 			      octave_idx_type& err, double& rcond,
-			      solve_singularity_handler sing_handler) const
+			      solve_singularity_handler sing_handler,
+			      bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -1411,23 +1446,26 @@
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
 	  octave_idx_type b_nc = b.cols ();
-	  rcond = 0.;
-
-	  // 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 += std::abs(data(i));
-	      if (atmp > anorm)
-		anorm = atmp;
+	  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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
 	    }
 
 	  if (typ == SparseType::Permuted_Upper)
 	    {
 	      retval.resize (nc, b_nc);
 	      octave_idx_type *perm = mattype.triangular_perm ();
-	      OCTAVE_LOCAL_BUFFER (Complex, work, nr);
+	      OCTAVE_LOCAL_BUFFER (Complex, work, nm);
 
 	      for (octave_idx_type j = 0; j < b_nc; j++)
 		{
@@ -1442,7 +1480,8 @@
 
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(kidx+1)-1) != k)
+			  if (ridx(cidx(kidx+1)-1) != k ||
+			      data(cidx(kidx+1)-1) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -1463,38 +1502,42 @@
 		    retval (perm[i], j) = work[i];
 		}
 
-	      // 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 (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++)
 		    {
-		      octave_idx_type iidx = perm[k];
-
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex 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 iidx = perm[k];
+
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type idx2 = ridx(i);
-			      work[idx2] = work[idx2] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  double atmp = 0;
-		  for (octave_idx_type i = 0; i < j+1; i++)
-		    {
-		      atmp += std::abs(work[i]);
-		      work[i] = 0.;
-		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
@@ -1513,7 +1556,8 @@
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k+1)-1) != k)
+			  if (ridx(cidx(k+1)-1) != k ||
+			      data(cidx(k+1)-1) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -1533,45 +1577,51 @@
 		    retval.xelem (i, j) = work[i];
 		}
 
-	      // 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 (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++)
 		    {
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex tmp = work[k] / data(cidx(k+1)-1);
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type iidx = ridx(i);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
-		      work[i] = 0.;
+		      double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
+			{
+			  atmp += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
-		}
-	    }
-
-	  rcond = 1. / ainvnorm / anorm;
+		  rcond = 1. / ainvnorm / anorm;
+		}
+	    }
 
 	triangular_error:
 	  if (err != 0)
 	    {
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -1585,7 +1635,10 @@
 	      err = -2;
 
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("matrix singular to machine precision, rcond = %g",
@@ -1602,7 +1655,8 @@
 SparseComplexMatrix
 SparseComplexMatrix::utsolve (SparseType &mattype, const SparseMatrix& b,
 			      octave_idx_type& err, double& rcond, 
-			      solve_singularity_handler sing_handler) const
+			      solve_singularity_handler sing_handler,
+			      bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -1625,20 +1679,23 @@
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
-	  rcond = 0.;
-
-	  // 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 += std::abs(data(i));
-	      if (atmp > anorm)
-		anorm = atmp;
+	  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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
 	    }
 
 	  octave_idx_type b_nc = b.cols ();
-	  octave_idx_type b_nz = b.nzmax ();
+	  octave_idx_type b_nz = b.nnz ();
 	  retval = SparseComplexMatrix (nc, b_nc, b_nz);
 	  retval.xcidx(0) = 0;
 	  octave_idx_type ii = 0;
@@ -1666,7 +1723,8 @@
 
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(kidx+1)-1) != k)
+			  if (ridx(cidx(kidx+1)-1) != k ||
+			      data(cidx(kidx+1)-1) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -1709,38 +1767,42 @@
 
 	      retval.maybe_compress ();
 
-	      // 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 (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++)
 		    {
-		      octave_idx_type iidx = perm[k];
-
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex 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 iidx = perm[k];
+
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type idx2 = ridx(i);
-			      work[idx2] = work[idx2] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  double atmp = 0;
-		  for (octave_idx_type i = 0; i < j+1; i++)
-		    {
-		      atmp += std::abs(work[i]);
-		      work[i] = 0.;
-		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
@@ -1758,7 +1820,8 @@
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k+1)-1) != k)
+			  if (ridx(cidx(k+1)-1) != k ||
+			      data(cidx(k+1)-1) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -1800,45 +1863,51 @@
 
 	      retval.maybe_compress ();
 
-	      // 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 (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++)
 		    {
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex tmp = work[k] / data(cidx(k+1)-1);
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type iidx = ridx(i);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
-		      work[i] = 0.;
+		      double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
+			{
+			  atmp += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
-		}
-	    }
-
-	  rcond = 1. / ainvnorm / anorm;
+		  rcond = 1. / ainvnorm / anorm;
+		}
+	    }
 
 	triangular_error:
 	  if (err != 0)
 	    {
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -1852,7 +1921,10 @@
 	      err = -2;
 
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("matrix singular to machine precision, rcond = %g",
@@ -1868,7 +1940,8 @@
 ComplexMatrix
 SparseComplexMatrix::utsolve (SparseType &mattype, const ComplexMatrix& b,
 			      octave_idx_type& err, double& rcond, 
-			      solve_singularity_handler sing_handler) const
+			      solve_singularity_handler sing_handler,
+			      bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -1892,16 +1965,19 @@
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
 	  octave_idx_type b_nc = b.cols ();
-	  rcond = 0.;
-
-	  // 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 += std::abs(data(i));
-	      if (atmp > anorm)
-		anorm = atmp;
+	  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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
 	    }
 
 	  if (typ == SparseType::Permuted_Upper)
@@ -1923,7 +1999,8 @@
 
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(kidx+1)-1) != k)
+			  if (ridx(cidx(kidx+1)-1) != k ||
+			      data(cidx(kidx+1)-1) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -1944,38 +2021,42 @@
 		    retval (perm[i], j) = work[i];
 		}
 
-	      // 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 (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++)
 		    {
-		      octave_idx_type iidx = perm[k];
-
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex 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 iidx = perm[k];
+
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type idx2 = ridx(i);
-			      work[idx2] = work[idx2] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  double atmp = 0;
-		  for (octave_idx_type i = 0; i < j+1; i++)
-		    {
-		      atmp += std::abs(work[i]);
-		      work[i] = 0.;
-		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
@@ -1994,7 +2075,8 @@
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k+1)-1) != k)
+			  if (ridx(cidx(k+1)-1) != k ||
+			      data(cidx(k+1)-1) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -2014,45 +2096,51 @@
 		    retval.xelem (i, j) = work[i];
 		}
 
-	      // 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 (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++)
 		    {
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex tmp = work[k] / data(cidx(k+1)-1);
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type iidx = ridx(i);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
-		      work[i] = 0.;
+		      double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
+			{
+			  atmp += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
-		}
-	    }
-
-	  rcond = 1. / ainvnorm / anorm;
+		  rcond = 1. / ainvnorm / anorm;
+		}
+	    }
 
 	triangular_error:
 	  if (err != 0)
 	    {
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -2066,7 +2154,10 @@
 	      err = -2;
 
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("matrix singular to machine precision, rcond = %g",
@@ -2083,7 +2174,8 @@
 SparseComplexMatrix
 SparseComplexMatrix::utsolve (SparseType &mattype, const SparseComplexMatrix& b,
 			      octave_idx_type& err, double& rcond, 
-			      solve_singularity_handler sing_handler) const
+			      solve_singularity_handler sing_handler,
+			      bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -2106,20 +2198,23 @@
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
-	  rcond = 0.;
-
-	  // 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 += std::abs(data(i));
-	      if (atmp > anorm)
-		anorm = atmp;
+	  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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
 	    }
 
 	  octave_idx_type b_nc = b.cols ();
-	  octave_idx_type b_nz = b.nzmax ();
+	  octave_idx_type b_nz = b.nnz ();
 	  retval = SparseComplexMatrix (nc, b_nc, b_nz);
 	  retval.xcidx(0) = 0;
 	  octave_idx_type ii = 0;
@@ -2147,7 +2242,8 @@
 
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(kidx+1)-1) != k)
+			  if (ridx(cidx(kidx+1)-1) != k ||
+			      data(cidx(kidx+1)-1) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -2190,38 +2286,42 @@
 
 	      retval.maybe_compress ();
 
-	      // 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 (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++)
 		    {
-		      octave_idx_type iidx = perm[k];
-
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex 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 iidx = perm[k];
+
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type idx2 = ridx(i);
-			      work[idx2] = work[idx2] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  double atmp = 0;
-		  for (octave_idx_type i = 0; i < j+1; i++)
-		    {
-		      atmp += std::abs(work[i]);
-		      work[i] = 0.;
-		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
@@ -2239,7 +2339,8 @@
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k+1)-1) != k)
+			  if (ridx(cidx(k+1)-1) != k ||
+			      data(cidx(k+1)-1) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -2281,45 +2382,51 @@
 
 	      retval.maybe_compress ();
 
-	      // 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 (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++)
 		    {
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex tmp = work[k] / data(cidx(k+1)-1);
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type iidx = ridx(i);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
-		      work[i] = 0.;
+		      double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
+			{
+			  atmp += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
-		}
-	    }
-
-	  rcond = 1. / ainvnorm / anorm;
+		  rcond = 1. / ainvnorm / anorm;
+		}
+	    }
 
 	triangular_error:
 	  if (err != 0)
 	    {
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -2333,7 +2440,10 @@
 	      err = -2;
 
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("matrix singular to machine precision, rcond = %g",
@@ -2350,7 +2460,8 @@
 ComplexMatrix
 SparseComplexMatrix::ltsolve (SparseType &mattype, const Matrix& b, 
 			      octave_idx_type& err, double& rcond, 
-			      solve_singularity_handler sing_handler) const
+			      solve_singularity_handler sing_handler,
+			      bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -2374,16 +2485,19 @@
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
 	  octave_idx_type b_nc = b.cols ();
-	  rcond = 0.;
-
-	  // 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 += std::abs(data(i));
-	      if (atmp > anorm)
-		anorm = atmp;
+	  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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
 	    }
 
 	  if (typ == SparseType::Permuted_Lower)
@@ -2413,7 +2527,7 @@
 				mini = i;
 			      }
 
-			  if (minr != k)
+			  if (minr != k || data (mini) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -2436,49 +2550,55 @@
 		    retval (i, j) = work[i];
 		}
 
-	      // 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 (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++)
 		    {
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = 0; k < nc; k++)
 			{
-			  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;
-			      }
-
-			  Complex tmp = work[k] / data(mini);
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
+			  if (work[k] != 0.)
 			    {
-			      if (i == mini)
-				continue;
-
-			      octave_idx_type iidx = perm[ridx(i)];
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
+				  }
+
+			      Complex 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 += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-
-		  double atmp = 0;
-		  for (octave_idx_type i = j; i < nc; i++)
-		    {
-		      atmp += std::abs(work[i]);
-		      work[i] = 0.;
-		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
@@ -2496,7 +2616,8 @@
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k)) != k)
+			  if (ridx(cidx(k)) != k ||
+			      data(cidx(k)) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -2515,46 +2636,51 @@
 		    retval.xelem (i, j) = work[i];
 		}
 
-	      // 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 (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++)
 		    {
-
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k < nc; k++)
 			{
-			  Complex tmp = work[k] / data(cidx(k));
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
+
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type iidx = ridx(i);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
-		      work[i] = 0.;
+		      double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
+			{
+			  atmp += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
-		}
-	    }
-
-	  rcond = 1. / ainvnorm / anorm;
-
+		  rcond = 1. / ainvnorm / anorm;
+		}
+	    }
 	triangular_error:
 	  if (err != 0)
 	    {
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -2568,7 +2694,10 @@
 	      err = -2;
 
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("matrix singular to machine precision, rcond = %g",
@@ -2585,7 +2714,8 @@
 SparseComplexMatrix
 SparseComplexMatrix::ltsolve (SparseType &mattype, const SparseMatrix& b, 
 			      octave_idx_type& err, double& rcond, 
-			      solve_singularity_handler sing_handler) const
+			      solve_singularity_handler sing_handler,
+			      bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -2609,20 +2739,23 @@
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
-	  rcond = 0.;
-
-	  // 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 += std::abs(data(i));
-	      if (atmp > anorm)
-		anorm = atmp;
+	  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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
 	    }
 
 	  octave_idx_type b_nc = b.cols ();
-	  octave_idx_type b_nz = b.nzmax ();
+	  octave_idx_type b_nz = b.nnz ();
 	  retval = SparseComplexMatrix (nc, b_nc, b_nz);
 	  retval.xcidx(0) = 0;
 	  octave_idx_type ii = 0;
@@ -2654,7 +2787,7 @@
 				mini = i;
 			      }
 
-			  if (minr != k)
+			  if (minr != k || data (mini) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -2699,49 +2832,55 @@
 
 	      retval.maybe_compress ();
 
-	      // 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 (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++)
 		    {
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = 0; k < nc; k++)
 			{
-			  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;
-			      }
-
-			  Complex tmp = work[k] / data(mini);
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
+			  if (work[k] != 0.)
 			    {
-			      if (i == mini)
-				continue;
-
-			      octave_idx_type iidx = perm[ridx(i)];
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
+				  }
+
+			      Complex 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 += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-
-		  double atmp = 0;
-		  for (octave_idx_type i = j; i < nc; i++)
-		    {
-		      atmp += std::abs(work[i]);
-		      work[i] = 0.;
-		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
@@ -2759,7 +2898,8 @@
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k)) != k)
+			  if (ridx(cidx(k)) != k ||
+			      data(cidx(k)) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -2801,47 +2941,52 @@
 
 	      retval.maybe_compress ();
 
-	      // 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 (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++)
 		    {
-
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k < nc; k++)
 			{
-			  Complex tmp = work[k] / data(cidx(k));
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
+
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type iidx = ridx(i);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
-		      work[i] = 0.;
+		      double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
+			{
+			  atmp += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
-		}
-
-	    }
-
-	  rcond = 1. / ainvnorm / anorm;
+		  rcond = 1. / ainvnorm / anorm;
+		}
+	    }
 
 	triangular_error:
 	  if (err != 0)
 	    {
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -2855,7 +3000,10 @@
 	      err = -2;
 
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("matrix singular to machine precision, rcond = %g",
@@ -2872,7 +3020,8 @@
 ComplexMatrix
 SparseComplexMatrix::ltsolve (SparseType &mattype, const ComplexMatrix& b,
 			      octave_idx_type& err, double& rcond,
-			      solve_singularity_handler sing_handler) const
+			      solve_singularity_handler sing_handler,
+			      bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -2896,16 +3045,19 @@
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
 	  octave_idx_type b_nc = b.cols ();
-	  rcond = 0.;
-
-	  // 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 += std::abs(data(i));
-	      if (atmp > anorm)
-		anorm = atmp;
+	  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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
 	    }
 
 	  if (typ == SparseType::Permuted_Lower)
@@ -2935,7 +3087,7 @@
 				mini = i;
 			      }
 
-			  if (minr != k)
+			  if (minr != k || data (mini) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -2958,49 +3110,55 @@
 		    retval (i, j) = work[i];
 		}
 
-	      // 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 (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++)
 		    {
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = 0; k < nc; k++)
 			{
-			  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;
-			      }
-
-			  Complex tmp = work[k] / data(mini);
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
+			  if (work[k] != 0.)
 			    {
-			      if (i == mini)
-				continue;
-
-			      octave_idx_type iidx = perm[ridx(i)];
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
+				  }
+
+			      Complex 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 += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-
-		  double atmp = 0;
-		  for (octave_idx_type i = j; i < nc; i++)
-		    {
-		      atmp += std::abs(work[i]);
-		      work[i] = 0.;
-		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
@@ -3020,7 +3178,8 @@
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k)) != k)
+			  if (ridx(cidx(k)) != k ||
+			      data(cidx(k)) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -3040,47 +3199,52 @@
 		    retval.xelem (i, j) = work[i];
 		}
 
-	      // 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 (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++)
 		    {
-
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k < nc; k++)
 			{
-			  Complex tmp = work[k] / data(cidx(k));
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
+
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type iidx = ridx(i);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
-		      work[i] = 0.;
+		      double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
+			{
+			  atmp += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
-		}
-
-	    }
-
-	  rcond = 1. / ainvnorm / anorm;
+		  rcond = 1. / ainvnorm / anorm;
+		}
+	    }
 
 	triangular_error:
 	  if (err != 0)
 	    {
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -3094,7 +3258,10 @@
 	      err = -2;
 
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("matrix singular to machine precision, rcond = %g",
@@ -3111,7 +3278,8 @@
 SparseComplexMatrix
 SparseComplexMatrix::ltsolve (SparseType &mattype, const SparseComplexMatrix& b,
 			      octave_idx_type& err, double& rcond, 
-			      solve_singularity_handler sing_handler) const
+			      solve_singularity_handler sing_handler,
+			      bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -3134,20 +3302,23 @@
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
-	  rcond = 0.;
-
-	  // 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 += std::abs(data(i));
-	      if (atmp > anorm)
-		anorm = atmp;
+	  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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
 	    }
 
 	  octave_idx_type b_nc = b.cols ();
-	  octave_idx_type b_nz = b.nzmax ();
+	  octave_idx_type b_nz = b.nnz ();
 	  retval = SparseComplexMatrix (nc, b_nc, b_nz);
 	  retval.xcidx(0) = 0;
 	  octave_idx_type ii = 0;
@@ -3179,7 +3350,7 @@
 				mini = i;
 			      }
 
-			  if (minr != k)
+			  if (minr != k || data (mini) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -3224,49 +3395,55 @@
 
 	      retval.maybe_compress ();
 
-	      // 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 (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++)
 		    {
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = 0; k < nc; k++)
 			{
-			  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;
-			      }
-
-			  Complex tmp = work[k] / data(mini);
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
+			  if (work[k] != 0.)
 			    {
-			      if (i == mini)
-				continue;
-
-			      octave_idx_type iidx = perm[ridx(i)];
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
+				  }
+
+			      Complex 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 += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-
-		  double atmp = 0;
-		  for (octave_idx_type i = j; i < nc; i++)
-		    {
-		      atmp += std::abs(work[i]);
-		      work[i] = 0.;
-		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
@@ -3284,7 +3461,8 @@
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k)) != k)
+			  if (ridx(cidx(k)) != k ||
+			      data(cidx(k)) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
@@ -3326,47 +3504,52 @@
 
 	      retval.maybe_compress ();
 
-	      // 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 (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++)
 		    {
-
-		      if (work[k] != 0.)
+		      work[j] = 1.;
+
+		      for (octave_idx_type k = j; k < nc; k++)
 			{
-			  Complex tmp = work[k] / data(cidx(k));
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
+
+			  if (work[k] != 0.)
 			    {
-			      octave_idx_type iidx = ridx(i);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      Complex 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 += std::abs(work[i]);
-		      work[i] = 0.;
+		      double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
+			{
+			  atmp += std::abs(work[i]);
+			  work[i] = 0.;
+			}
+		      if (atmp > ainvnorm)
+			ainvnorm = atmp;
 		    }
-		  if (atmp > ainvnorm)
-		    ainvnorm = atmp;
-		}
-
-	    }
-
-	  rcond = 1. / ainvnorm / anorm;
+		  rcond = 1. / ainvnorm / anorm;
+		}
+	    }
 
 	triangular_error:
 	  if (err != 0)
 	    {
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -3380,7 +3563,10 @@
 	      err = -2;
 
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("matrix singular to machine precision, rcond = %g",
@@ -3395,9 +3581,10 @@
 }
 
 ComplexMatrix
-SparseComplexMatrix::trisolve (SparseType &mattype, const Matrix& b, octave_idx_type& err, 
-			       double& rcond,
-			       solve_singularity_handler sing_handler) const
+SparseComplexMatrix::trisolve (SparseType &mattype, const Matrix& b,
+			       octave_idx_type& err, double& rcond,
+			       solve_singularity_handler sing_handler,
+			       bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -3408,6 +3595,9 @@
   if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
+  else if (calc_cond)
+    (*current_liboctave_error_handler) 
+      ("calculation of condition number not implemented");
   else
     {
       // Print spparms("spumoni") info if requested
@@ -3526,7 +3716,10 @@
 	      err = -2;
 
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("matrix singular to machine precision");
@@ -3544,8 +3737,9 @@
 
 SparseComplexMatrix
 SparseComplexMatrix::trisolve (SparseType &mattype, const SparseMatrix& b,
-		     octave_idx_type& err, double& rcond, 
-		     solve_singularity_handler sing_handler) const
+			       octave_idx_type& err, double& rcond, 
+			       solve_singularity_handler sing_handler,
+			       bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -3556,6 +3750,9 @@
   if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
+  else if (calc_cond)
+    (*current_liboctave_error_handler) 
+      ("calculation of condition number not implemented");
   else
     {
       // Print spparms("spumoni") info if requested
@@ -3614,13 +3811,16 @@
 	      ("unrecoverable error in zgttrf");
 	  else
 	    {
-	      rcond = 0.0;
 	      if (err != 0) 
 		{
 		  err = -2;
+		  rcond = 0.0;
 
 		  if (sing_handler)
-		    sing_handler (rcond);
+		    {
+		      sing_handler (rcond);
+		      mattype.mark_as_rectangular ();
+		    }
 		  else
 		    (*current_liboctave_error_handler)
 		      ("matrix singular to machine precision");
@@ -3629,11 +3829,12 @@
 	      else 
 		{
 		  char job = 'N';
-		  volatile octave_idx_type x_nz = b.nzmax ();
+		  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;
+		  rcond = 1.0;
 
 		  OCTAVE_LOCAL_BUFFER (Complex, work, nr);
 
@@ -3695,7 +3896,8 @@
 ComplexMatrix
 SparseComplexMatrix::trisolve (SparseType &mattype, const ComplexMatrix& b,
 			       octave_idx_type& err, double& rcond, 
-			       solve_singularity_handler sing_handler) const
+			       solve_singularity_handler sing_handler,
+			       bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -3706,6 +3908,9 @@
   if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
+  else if (calc_cond)
+    (*current_liboctave_error_handler) 
+      ("calculation of condition number not implemented");
   else
     {
       // Print spparms("spumoni") info if requested
@@ -3834,7 +4039,10 @@
 	      err = -2;
 		      
 	      if (sing_handler)
-		sing_handler (rcond);
+		{
+		  sing_handler (rcond);
+		  mattype.mark_as_rectangular ();
+		}
 	      else
 		(*current_liboctave_error_handler)
 		  ("matrix singular to machine precision");
@@ -3849,8 +4057,10 @@
 
 SparseComplexMatrix
 SparseComplexMatrix::trisolve (SparseType &mattype, 
-		     const SparseComplexMatrix& b, octave_idx_type& err, double& rcond, 
-		     solve_singularity_handler sing_handler) const
+			       const SparseComplexMatrix& b, 
+			       octave_idx_type& err, double& rcond, 
+			       solve_singularity_handler sing_handler,
+			       bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -3861,6 +4071,9 @@
   if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
+  else if (calc_cond)
+    (*current_liboctave_error_handler) 
+      ("calculation of condition number not implemented");
   else
     {
       // Print spparms("spumoni") info if requested
@@ -3919,13 +4132,16 @@
 	      ("unrecoverable error in zgttrf");
 	  else
 	    {
-	      rcond = 0.0;
 	      if (err != 0) 
 		{
+		  rcond = 0.0;
 		  err = -2;
 
 		  if (sing_handler)
-		    sing_handler (rcond);
+		    {
+		      sing_handler (rcond);
+		      mattype.mark_as_rectangular ();
+		    }
 		  else
 		    (*current_liboctave_error_handler)
 		      ("matrix singular to machine precision");
@@ -3940,7 +4156,7 @@
 
 		  // 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.nzmax ();
+		  volatile octave_idx_type x_nz = b.nnz ();
 		  volatile octave_idx_type ii = 0;
 		  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
 
@@ -4010,9 +4226,10 @@
 }
 
 ComplexMatrix
-SparseComplexMatrix::bsolve (SparseType &mattype, const Matrix& b, octave_idx_type& err, 
-			     double& rcond,
-			     solve_singularity_handler sing_handler) const
+SparseComplexMatrix::bsolve (SparseType &mattype, const Matrix& b,
+			     octave_idx_type& err, double& rcond,
+			     solve_singularity_handler sing_handler,
+			     bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -4054,7 +4271,9 @@
 	      }
 
 	  // Calculate the norm of the matrix, for later use.
-	  // double anorm = m_band.abs().sum().row(0).max();
+	  double anorm;
+	  if (calc_cond)
+	    anorm = m_band.abs().sum().row(0).max();
 
 	  char job = 'L';
 	  F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
@@ -4066,9 +4285,9 @@
 	      ("unrecoverable error in zpbtrf");
 	  else
 	    {
-	      rcond = 0.0;
 	      if (err != 0) 
 		{
+		  rcond = 0.0;
 		  // Matrix is not positive definite!! Fall through to
 		  // unsymmetric banded solver.
 		  mattype.mark_as_unsymmetric ();
@@ -4077,68 +4296,69 @@
 		} 
 	      else 
 		{
-		  // Unfortunately, the time to calculate the condition
-		  // number is dominant for narrow banded matrices and
-		  // so we rely on the "err" flag from xPBTRF to flag
-		  // singularity. The commented code below is left here
-		  // for reference
-
-		  //Array<double> z (3 * nr);
-		  //Complex *pz = z.fortran_vec ();
-		  //Array<octave_idx_type> iz (nr);
-		  //octave_idx_type *piz = iz.fortran_vec ();
-		  //
-		  //F77_XFCN (zpbcon, ZGBCON, 
-		  //	(F77_CONST_CHAR_ARG2 (&job, 1),
-		  //	 nr, n_lower, tmp_data, ldm,
-		  //	 anorm, rcond, pz, piz, err
-		  //	 F77_CHAR_ARG_LEN (1)));
-		  //
-		  //
-		  //if (f77_exception_encountered)
-		  //	(*current_liboctave_error_handler) 
-		  //	  ("unrecoverable error in zpbcon");
-		  //
-		  //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);
-		  //    else
-		  //      (*current_liboctave_error_handler)
-		  //        ("matrix singular to machine precision, rcond = %g",
-		  //         rcond);
-		  //  }
-		  //else
-		  //    REST OF CODE, EXCEPT rcond=1
-
-		  rcond = 1.;
-		  retval = ComplexMatrix (b);
-		  Complex *result = retval.fortran_vec ();
-
-		  octave_idx_type b_nc = b.cols ();
-
-		  F77_XFCN (zpbtrs, ZPBTRS, 
-			    (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nr, n_lower, b_nc, tmp_data,
-			     ldm, result, b.rows(), err
-			     F77_CHAR_ARG_LEN (1)));
+		  if (calc_cond)
+		    {
+		      Array<Complex> z (2 * nr);
+		      Complex *pz = z.fortran_vec ();
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
+
+		      F77_XFCN (zpbcon, ZPBCON, 
+		      	(F77_CONST_CHAR_ARG2 (&job, 1),
+		      	 nr, n_lower, tmp_data, ldm,
+		      	 anorm, rcond, pz, piz, err
+		      	 F77_CHAR_ARG_LEN (1)));
+
+		      if (f77_exception_encountered)
+		      	(*current_liboctave_error_handler) 
+		      	  ("unrecoverable error in zpbcon");
+
+		      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.0;
+
+		  if (err == 0)
+		    {
+		      retval = ComplexMatrix (b);
+		      Complex *result = retval.fortran_vec ();
+
+		      octave_idx_type b_nc = b.cols ();
+
+		      F77_XFCN (zpbtrs, ZPBTRS, 
+				(F77_CONST_CHAR_ARG2 (&job, 1),
+				 nr, n_lower, b_nc, tmp_data,
+				 ldm, result, b.rows(), err
+				 F77_CHAR_ARG_LEN (1)));
 		    
-		  if (f77_exception_encountered)
-		    (*current_liboctave_error_handler)
-		      ("unrecoverable error in zpbtrs");
-
-		  if (err != 0)
-		    {
-		      (*current_liboctave_error_handler) 
-			("SparseMatrix::solve solve failed");
-		      err = -1;
+		      if (f77_exception_encountered)
+			(*current_liboctave_error_handler)
+			  ("unrecoverable error in zpbtrs");
+
+		      if (err != 0)
+			{
+			  (*current_liboctave_error_handler) 
+			    ("SparseMatrix::solve solve failed");
+			  err = -1;
+			}
 		    }
 		}
 	    }
@@ -4167,6 +4387,20 @@
 	    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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
+	    }
+
 	  Array<octave_idx_type> ipvt (nr);
 	  octave_idx_type *pipvt = ipvt.fortran_vec ();
 
@@ -4180,73 +4414,81 @@
 	    {
 	      // Throw-away extra info LAPACK gives so as to not 
 	      // change output.
-	      rcond = 0.0;
 	      if (err != 0) 
 		{
+		  rcond = 0.0;
 		  err = -2;
 
 		  if (sing_handler)
-		    sing_handler (rcond);
+		    {
+		      sing_handler (rcond);
+		      mattype.mark_as_rectangular ();
+		    }
 		  else
 		    (*current_liboctave_error_handler)
 		      ("matrix singular to machine precision");
-
 		} 
 	      else 
 		{
-		  char job = '1';
-
-		  // Unfortunately, the time to calculate the condition
-		  // number is dominant for narrow banded matrices and
-		  // so we rely on the "err" flag from xPBTRF to flag
-		  // singularity. The commented code below is left here
-		  // for reference
-
-		  //F77_XFCN (zgbcon, ZGBCON, 
-		  //	(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 (f77_exception_encountered)
-		  //  (*current_liboctave_error_handler) 
-		  //    ("unrecoverable error in zgbcon");
-		  //
-		  // 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);
-		  //    else
-		  //      (*current_liboctave_error_handler)
-		  //        ("matrix singular to machine precision, rcond = %g",
-		  //         rcond);
-		  //  }
-		  //else
-		  //  REST OF CODE, EXCEPT rcond=1
-
-		  rcond = 1.;
-		  retval = ComplexMatrix (b);
-		  Complex *result = retval.fortran_vec ();
-
-		  octave_idx_type b_nc = b.cols ();
-
-		  job = 'N';
-		  F77_XFCN (zgbtrs, ZGBTRS, 
-			    (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)));
+		  if (calc_cond)
+		    {
+		      char job = '1';
+		      Array<Complex> z (2 * nr);
+		      Complex *pz = z.fortran_vec ();
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
+
+		      F77_XFCN (zgbcon, ZGBCON, 
+		      	(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 (f77_exception_encountered)
+		        (*current_liboctave_error_handler) 
+		          ("unrecoverable error in zgbcon");
+
+		       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 = ComplexMatrix (b);
+		      Complex *result = retval.fortran_vec ();
+
+		      octave_idx_type b_nc = b.cols ();
+
+		      char job = 'N';
+		      F77_XFCN (zgbtrs, ZGBTRS, 
+				(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)));
 		    
-		  if (f77_exception_encountered)
-		    (*current_liboctave_error_handler)
-		      ("unrecoverable error in zgbtrs");
+		      if (f77_exception_encountered)
+			(*current_liboctave_error_handler)
+			  ("unrecoverable error in zgbtrs");
+		    }
 		}
 	    }
 	}
@@ -4260,7 +4502,8 @@
 SparseComplexMatrix
 SparseComplexMatrix::bsolve (SparseType &mattype, const SparseMatrix& b,
 			     octave_idx_type& err, double& rcond, 
-			     solve_singularity_handler sing_handler) const
+			     solve_singularity_handler sing_handler,
+			     bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -4302,6 +4545,11 @@
 		  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 (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
 				     nr, n_lower, tmp_data, ldm, err
@@ -4312,75 +4560,118 @@
 	      ("unrecoverable error in zpbtrf");
 	  else
 	    {
-	      rcond = 0.0;
 	      if (err != 0) 
 		{
+		  rcond = 0.0;
 		  mattype.mark_as_unsymmetric ();
 		  typ = SparseType::Banded;
 		  err = 0;
 		} 
 	      else 
 		{
-		  rcond = 1.;
-		  octave_idx_type b_nr = b.rows ();
-		  octave_idx_type b_nc = b.cols ();
-		  OCTAVE_LOCAL_BUFFER (Complex, 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.nzmax ();
-		  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++)
+		  if (calc_cond)
+		    {
+		      Array<Complex> z (2 * nr);
+		      Complex *pz = z.fortran_vec ();
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
+
+		      F77_XFCN (zpbcon, ZPBCON, 
+		      	(F77_CONST_CHAR_ARG2 (&job, 1),
+		      	 nr, n_lower, tmp_data, ldm,
+		      	 anorm, rcond, pz, piz, err
+		      	 F77_CHAR_ARG_LEN (1)));
+
+		      if (f77_exception_encountered)
+		      	(*current_liboctave_error_handler) 
+		      	  ("unrecoverable error in zpbcon");
+
+		      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.0;
+
+		  if (err == 0)
 		    {
-		      for (octave_idx_type i = 0; i < b_nr; i++)
-			Bx[i] = b.elem (i, j);
-
-		      F77_XFCN (zpbtrs, ZPBTRS, 
-				(F77_CONST_CHAR_ARG2 (&job, 1),
-				 nr, n_lower, 1, tmp_data,
-				 ldm, Bx, b_nr, err
-				 F77_CHAR_ARG_LEN (1)));
+		      octave_idx_type b_nr = b.rows ();
+		      octave_idx_type b_nc = b.cols ();
+		      OCTAVE_LOCAL_BUFFER (Complex, 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 = 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++)
+			    Bx[i] = b.elem (i, j);
+
+			  F77_XFCN (zpbtrs, ZPBTRS, 
+				    (F77_CONST_CHAR_ARG2 (&job, 1),
+				     nr, n_lower, 1, tmp_data,
+				     ldm, Bx, b_nr, err
+				     F77_CHAR_ARG_LEN (1)));
 		    
-		      if (f77_exception_encountered)
-			{
-			  (*current_liboctave_error_handler)
-			    ("unrecoverable error in dpbtrs");
-			  err = -1;
-			  break;
+			  if (f77_exception_encountered)
+			    {
+			      (*current_liboctave_error_handler)
+				("unrecoverable error in dpbtrs");
+			      err = -1;
+			      break;
+			    }
+
+			  if (err != 0)
+			    {
+			      (*current_liboctave_error_handler) 
+				("SparseComplexMatrix::solve solve failed");
+			      err = -1;
+			      break;
+			    }
+
+			  for (octave_idx_type i = 0; i < b_nr; i++)
+			    {
+			      Complex 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;
 			}
 
-		      if (err != 0)
-			{
-			  (*current_liboctave_error_handler) 
-			    ("SparseComplexMatrix::solve solve failed");
-			  err = -1;
-			  break;
-			}
-
-		      for (octave_idx_type i = 0; i < b_nr; i++)
-			{
-			  Complex 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 ();
 		    }
-
-		  retval.maybe_compress ();
 		}
 	    }
 	}
@@ -4408,6 +4699,20 @@
 	    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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
+	    }
+
 	  Array<octave_idx_type> ipvt (nr);
 	  octave_idx_type *pipvt = ipvt.fortran_vec ();
 
@@ -4419,13 +4724,16 @@
 	      ("unrecoverable error in zgbtrf");
 	  else
 	    {
-	      rcond = 0.0;
 	      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");
@@ -4433,60 +4741,105 @@
 		} 
 	      else 
 		{
-		  char job = 'N';
-		  volatile octave_idx_type x_nz = b.nzmax ();
-		  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 (Complex, work, nr);
-
-		  for (volatile octave_idx_type j = 0; j < b_nc; j++)
+		  if (calc_cond)
 		    {
-		      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 (zgbtrs, ZGBTRS, 
-				(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)));
-		    
+		      char job = '1';
+		      Array<Complex> z (2 * nr);
+		      Complex *pz = z.fortran_vec ();
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
+
+		      F77_XFCN (zgbcon, ZGBCON, 
+		      	(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 (f77_exception_encountered)
+		        (*current_liboctave_error_handler) 
+		          ("unrecoverable error in zgbcon");
+
+		       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 (Complex, work, nr);
+
+		      for (volatile octave_idx_type j = 0; j < b_nc; j++)
 			{
-			  (*current_liboctave_error_handler)
-			    ("unrecoverable error in zgbtrs");
-			  break;
+			  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 (zgbtrs, ZGBTRS, 
+				    (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)));
+		    
+			  if (f77_exception_encountered)
+			    {
+			      (*current_liboctave_error_handler)
+				("unrecoverable error in zgbtrs");
+			      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 (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;
 			}
 
-		      // 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 ();
 		    }
-
-		  retval.maybe_compress ();
 		}
 	    }
 	}
@@ -4500,7 +4853,8 @@
 ComplexMatrix
 SparseComplexMatrix::bsolve (SparseType &mattype, const ComplexMatrix& b, 
 			     octave_idx_type& err, double& rcond, 
-			     solve_singularity_handler sing_handler) const
+			     solve_singularity_handler sing_handler,
+			     bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -4542,6 +4896,11 @@
 		  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 (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
 				     nr, n_lower, tmp_data, ldm, err
@@ -4552,41 +4911,83 @@
 	      ("unrecoverable error in zpbtrf");
 	  else
 	    {
-	      rcond = 0.0;
 	      if (err != 0) 
 		{
 		  // Matrix is not positive definite!! Fall through to
 		  // unsymmetric banded solver.
+		  rcond = 0.0;
 		  mattype.mark_as_unsymmetric ();
 		  typ = SparseType::Banded;
 		  err = 0;
 		} 
 	      else 
 		{
-		  rcond = 1.;
-		  octave_idx_type b_nr = b.rows ();
-		  octave_idx_type b_nc = b.cols ();
-		  retval = ComplexMatrix (b);
-		  Complex *result = retval.fortran_vec ();
-
-		  F77_XFCN (zpbtrs, ZPBTRS, 
-			    (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nr, n_lower, b_nc, tmp_data,
-			     ldm, result, b_nr, err
-			     F77_CHAR_ARG_LEN (1)));
+		  if (calc_cond)
+		    {
+		      Array<Complex> z (2 * nr);
+		      Complex *pz = z.fortran_vec ();
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
+
+		      F77_XFCN (zpbcon, ZPBCON, 
+		      	(F77_CONST_CHAR_ARG2 (&job, 1),
+		      	 nr, n_lower, tmp_data, ldm,
+		      	 anorm, rcond, pz, piz, err
+		      	 F77_CHAR_ARG_LEN (1)));
+
+		      if (f77_exception_encountered)
+		      	(*current_liboctave_error_handler) 
+		      	  ("unrecoverable error in zpbcon");
+
+		      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.0;
+
+		  if (err == 0)
+		    {
+		      octave_idx_type b_nr = b.rows ();
+		      octave_idx_type b_nc = b.cols ();
+		      retval = ComplexMatrix (b);
+		      Complex *result = retval.fortran_vec ();
+
+		      F77_XFCN (zpbtrs, ZPBTRS, 
+				(F77_CONST_CHAR_ARG2 (&job, 1),
+				 nr, n_lower, b_nc, tmp_data,
+				 ldm, result, b_nr, err
+				 F77_CHAR_ARG_LEN (1)));
 		    
-		  if (f77_exception_encountered)
-		    {
-		      (*current_liboctave_error_handler)
-			("unrecoverable error in zpbtrs");
-		      err = -1;
-		    }
-
-		  if (err != 0)
-		    {
-		      (*current_liboctave_error_handler) 
-			("SparseComplexMatrix::solve solve failed");
-		      err = -1;
+		      if (f77_exception_encountered)
+			{
+			  (*current_liboctave_error_handler)
+			    ("unrecoverable error in zpbtrs");
+			  err = -1;
+			}
+
+		      if (err != 0)
+			{
+			  (*current_liboctave_error_handler) 
+			    ("SparseComplexMatrix::solve solve failed");
+			  err = -1;
+			}
 		    }
 		}
 	    }
@@ -4615,6 +5016,20 @@
 	    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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
+	    }
+
 	  Array<octave_idx_type> ipvt (nr);
 	  octave_idx_type *pipvt = ipvt.fortran_vec ();
 
@@ -4626,35 +5041,81 @@
 	      ("unrecoverable error in zgbtrf");
 	  else
 	    {
-	      rcond = 0.0;
 	      if (err != 0) 
 		{
 		  err = -2;
+		  rcond = 0.0;
 
 		  if (sing_handler)
-		    sing_handler (rcond);
+		    {
+		      sing_handler (rcond);
+		      mattype.mark_as_rectangular ();
+		    }
 		  else
 		    (*current_liboctave_error_handler)
 		      ("matrix singular to machine precision");
-
 		} 
 	      else 
 		{
-		  char job = 'N';
-		  octave_idx_type b_nc = b.cols ();
-		  retval = ComplexMatrix (b);
-		  Complex *result = retval.fortran_vec ();
-
-		  F77_XFCN (zgbtrs, ZGBTRS, 
-			    (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)));
+		  if (calc_cond)
+		    {
+		      char job = '1';
+		      Array<Complex> z (2 * nr);
+		      Complex *pz = z.fortran_vec ();
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
+
+		      F77_XFCN (zgbcon, ZGBCON, 
+		      	(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 (f77_exception_encountered)
+		        (*current_liboctave_error_handler) 
+		          ("unrecoverable error in zgbcon");
+
+		       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 = ComplexMatrix (b);
+		      Complex *result = retval.fortran_vec ();
+
+		      F77_XFCN (zgbtrs, ZGBTRS, 
+				(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)));
 		    
-		  if (f77_exception_encountered)
-		    {
-		      (*current_liboctave_error_handler)
-			("unrecoverable error in dgbtrs");
+		      if (f77_exception_encountered)
+			{
+			  (*current_liboctave_error_handler)
+			    ("unrecoverable error in dgbtrs");
+			}
 		    }
 		}
 	    }
@@ -4668,8 +5129,9 @@
 
 SparseComplexMatrix
 SparseComplexMatrix::bsolve (SparseType &mattype, const SparseComplexMatrix& b,
-		     octave_idx_type& err, double& rcond, 
-		     solve_singularity_handler sing_handler) const
+			     octave_idx_type& err, double& rcond, 
+			     solve_singularity_handler sing_handler,
+			     bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -4711,6 +5173,11 @@
 		  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 (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
 				     nr, n_lower, tmp_data, ldm, err
@@ -4721,7 +5188,6 @@
 	      ("unrecoverable error in zpbtrf");
 	  else
 	    {
-	      rcond = 0.0;
 	      if (err != 0) 
 		{
 		  // Matrix is not positive definite!! Fall through to
@@ -4729,77 +5195,119 @@
 		  mattype.mark_as_unsymmetric ();
 		  typ = SparseType::Banded;
 
+		  rcond = 0.0;
 		  err = 0;
 		} 
 	      else 
 		{
-		  rcond = 1.;
-		  octave_idx_type b_nr = b.rows ();
-		  octave_idx_type b_nc = b.cols ();
-		  OCTAVE_LOCAL_BUFFER (Complex, 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.nzmax ();
-		  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++)
+		  if (calc_cond)
 		    {
-
-		      for (octave_idx_type i = 0; i < b_nr; i++)
-			Bx[i] = b (i,j);
-
-		      F77_XFCN (zpbtrs, ZPBTRS, 
-				(F77_CONST_CHAR_ARG2 (&job, 1),
-				 nr, n_lower, 1, tmp_data,
-				 ldm, Bx, b_nr, err
-				 F77_CHAR_ARG_LEN (1)));
-		    
+		      Array<Complex> z (2 * nr);
+		      Complex *pz = z.fortran_vec ();
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
+
+		      F77_XFCN (zpbcon, ZPBCON, 
+		      	(F77_CONST_CHAR_ARG2 (&job, 1),
+		      	 nr, n_lower, tmp_data, ldm,
+		      	 anorm, rcond, pz, piz, err
+		      	 F77_CHAR_ARG_LEN (1)));
+
 		      if (f77_exception_encountered)
-			{
-			  (*current_liboctave_error_handler)
-			    ("unrecoverable error in zpbtrs");
-			  err = -1;
-			  break;
-			}
-
-		      if (err != 0)
+		      	(*current_liboctave_error_handler) 
+		      	  ("unrecoverable error in zpbcon");
+
+		      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.0;
+
+		  if (err == 0)
+		    {
+		      octave_idx_type b_nr = b.rows ();
+		      octave_idx_type b_nc = b.cols ();
+		      OCTAVE_LOCAL_BUFFER (Complex, 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 = SparseComplexMatrix (b_nr, b_nc, x_nz);
+
+		      retval.xcidx(0) = 0;
+		      for (volatile octave_idx_type j = 0; j < b_nc; j++)
 			{
-			  (*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.)
-			  new_nnz++;
+
+			  for (octave_idx_type i = 0; i < b_nr; i++)
+			    Bx[i] = b (i,j);
+
+			  F77_XFCN (zpbtrs, ZPBTRS, 
+				    (F77_CONST_CHAR_ARG2 (&job, 1),
+				     nr, n_lower, 1, tmp_data,
+				     ldm, Bx, b_nr, err
+				     F77_CHAR_ARG_LEN (1)));
+		    
+			  if (f77_exception_encountered)
+			    {
+			      (*current_liboctave_error_handler)
+				("unrecoverable error in zpbtrs");
+			      err = -1;
+			      break;
+			    }
+
+			  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.)
+			      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;
-			}
+			  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.)
-			  {
-			    retval.xridx(ii) = i;
-			    retval.xdata(ii++) = Bx[i];
-			  }
-
-		      retval.xcidx(j+1) = ii;
+			  for (octave_idx_type i = 0; i < nr; i++)
+			    if (Bx[i] != 0.)
+			      {
+				retval.xridx(ii) = i;
+				retval.xdata(ii++) = Bx[i];
+			      }
+
+			  retval.xcidx(j+1) = ii;
+			}
+
+		      retval.maybe_compress ();
 		    }
-
-		  retval.maybe_compress ();
 		}
 	    }
 	}
@@ -4827,6 +5335,20 @@
 	    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 += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
+	    }
+
 	  Array<octave_idx_type> ipvt (nr);
 	  octave_idx_type *pipvt = ipvt.fortran_vec ();
 
@@ -4838,13 +5360,16 @@
 	      ("unrecoverable error in xgbtrf");
 	  else
 	    {
-	      rcond = 0.0;
 	      if (err != 0) 
 		{
 		  err = -2;
+		  rcond = 0.0;
 
 		  if (sing_handler)
-		    sing_handler (rcond);
+		    {
+		      sing_handler (rcond);
+		      mattype.mark_as_rectangular ();
+		    }
 		  else
 		    (*current_liboctave_error_handler)
 		      ("matrix singular to machine precision");
@@ -4852,61 +5377,106 @@
 		}
 	      else 
 		{
-		  char job = 'N';
-		  volatile octave_idx_type x_nz = b.nzmax ();
-		  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 (Complex, Bx, nr);
-
-		  for (volatile octave_idx_type j = 0; j < b_nc; j++)
+		  if (calc_cond)
 		    {
-		      for (octave_idx_type i = 0; i < nr; i++)
-			Bx[i] = 0.;
-
-		      for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
-			Bx[b.ridx(i)] = b.data(i);
-
-		      F77_XFCN (zgbtrs, ZGBTRS, 
-				(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)));
-		    
+		      char job = '1';
+		      Array<Complex> z (2 * nr);
+		      Complex *pz = z.fortran_vec ();
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
+
+		      F77_XFCN (zgbcon, ZGBCON, 
+		      	(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 (f77_exception_encountered)
+		        (*current_liboctave_error_handler) 
+		          ("unrecoverable error in zgbcon");
+
+		       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 (Complex, Bx, nr);
+
+		      for (volatile octave_idx_type j = 0; j < b_nc; j++)
 			{
-			  (*current_liboctave_error_handler)
-			    ("unrecoverable error in dgbtrs");
-			  break;
+			  for (octave_idx_type i = 0; i < nr; i++)
+			    Bx[i] = 0.;
+
+			  for (octave_idx_type i = b.cidx(j); 
+			       i < b.cidx(j+1); i++)
+			    Bx[b.ridx(i)] = b.data(i);
+
+			  F77_XFCN (zgbtrs, ZGBTRS, 
+				    (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)));
+		    
+			  if (f77_exception_encountered)
+			    {
+			      (*current_liboctave_error_handler)
+				("unrecoverable error in dgbtrs");
+			      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.)
+			      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.)
+			      {
+				retval.xridx(ii) = i;
+				retval.xdata(ii++) = Bx[i]; 
+			      }
+			  retval.xcidx(j+1) = ii;
 			}
 
-		      // 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.)
-			  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.)
-			  {
-			    retval.xridx(ii) = i;
-			    retval.xdata(ii++) = Bx[i]; 
-			  }
-		      retval.xcidx(j+1) = ii;
+		      retval.maybe_compress ();
 		    }
-
-		  retval.maybe_compress ();
 		}
 	    }
 	}
@@ -4918,9 +5488,10 @@
 }
 
 void *
-SparseComplexMatrix::factorize (octave_idx_type& err, double &rcond, Matrix &Control, 
-				Matrix &Info,
-				solve_singularity_handler sing_handler) const
+SparseComplexMatrix::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;
@@ -4985,7 +5556,10 @@
 				   Symbolic, &Numeric, control, info) ;
       UMFPACK_ZNAME (free_symbolic) (&Symbolic) ;
 
-      rcond = Info (UMFPACK_RCOND);
+      if (calc_cond)
+	rcond = Info (UMFPACK_RCOND);
+      else
+	rcond = 1.;
       volatile double rcond_plus_one = rcond + 1.0;
 
       if (status == UMFPACK_WARNING_singular_matrix || 
@@ -5029,9 +5603,10 @@
 }
 
 ComplexMatrix
-SparseComplexMatrix::fsolve (SparseType &mattype, const Matrix& b, octave_idx_type& err, 
-			     double& rcond,
-			     solve_singularity_handler sing_handler) const
+SparseComplexMatrix::fsolve (SparseType &mattype, const Matrix& b,
+			     octave_idx_type& err, double& rcond,
+			     solve_singularity_handler sing_handler,
+			     bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -5139,7 +5714,10 @@
 	  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 	  L = CHOLMOD_NAME(analyze) (A, cm);
 	  CHOLMOD_NAME(factorize) (A, L, cm);
-	  rcond = CHOLMOD_NAME(rcond)(L, cm);
+	  if (calc_cond)
+	    rcond = CHOLMOD_NAME(rcond)(L, cm);
+	  else
+	    rcond = 1.;
 	  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 
 	  if (rcond == 0.0)
@@ -5157,7 +5735,10 @@
 		  err = -2;
 
 		  if (sing_handler)
-		    sing_handler (rcond);
+		    {
+		      sing_handler (rcond);
+		      mattype.mark_as_rectangular ();
+		    }
 		  else
 		    (*current_liboctave_error_handler)
 		      ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -5200,7 +5781,7 @@
 #ifdef HAVE_UMFPACK
 	  Matrix Control, Info;
 	  void *Numeric = factorize (err, rcond, Control, Info, 
-				     sing_handler);
+				     sing_handler, calc_cond);
 
 	  if (err == 0)
 	    {
@@ -5264,6 +5845,9 @@
 
 	      UMFPACK_ZNAME (free_numeric) (&Numeric);
 	    }
+	  else
+	    mattype.mark_as_rectangular ();
+
 #else
 	  (*current_liboctave_error_handler) ("UMFPACK not installed");
 #endif
@@ -5278,7 +5862,8 @@
 SparseComplexMatrix
 SparseComplexMatrix::fsolve (SparseType &mattype, const SparseMatrix& b, 
 			     octave_idx_type& err, double& rcond,
-			     solve_singularity_handler sing_handler) const
+			     solve_singularity_handler sing_handler,
+			     bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -5397,7 +5982,10 @@
 	  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 	  L = CHOLMOD_NAME(analyze) (A, cm);
 	  CHOLMOD_NAME(factorize) (A, L, cm);
-	  rcond = CHOLMOD_NAME(rcond)(L, cm);
+	  if (calc_cond)
+	    rcond = CHOLMOD_NAME(rcond)(L, cm);
+	  else
+	    rcond = 1.;
 	  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 
 	  if (rcond == 0.0)
@@ -5415,7 +6003,10 @@
 		  err = -2;
 
 		  if (sing_handler)
-		    sing_handler (rcond);
+		    {
+		      sing_handler (rcond);
+		      mattype.mark_as_rectangular ();
+		    }
 		  else
 		    (*current_liboctave_error_handler)
 		      ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -5463,7 +6054,8 @@
 	{
 #ifdef HAVE_UMFPACK
 	  Matrix Control, Info;
-	  void *Numeric = factorize (err, rcond, Control, Info, sing_handler);
+	  void *Numeric = factorize (err, rcond, Control, Info, 
+				     sing_handler, calc_cond);
 
 	  if (err == 0)
 	    {
@@ -5487,7 +6079,7 @@
 
 	      // Take a first guess that the number of non-zero terms
 	      // will be as many as in b
-	      octave_idx_type x_nz = b.nzmax ();
+	      octave_idx_type x_nz = b.nnz ();
 	      octave_idx_type ii = 0;
 	      retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
 
@@ -5557,6 +6149,9 @@
 
 	      UMFPACK_ZNAME (free_numeric) (&Numeric);
 	    }
+	  else
+	    mattype.mark_as_rectangular ();
+
 #else
 	  (*current_liboctave_error_handler) ("UMFPACK not installed");
 #endif
@@ -5571,7 +6166,8 @@
 ComplexMatrix
 SparseComplexMatrix::fsolve (SparseType &mattype, const ComplexMatrix& b, 
 			     octave_idx_type& err, double& rcond,
-			     solve_singularity_handler sing_handler) const
+			     solve_singularity_handler sing_handler,
+			     bool calc_cond) const
 {
   ComplexMatrix retval;
 
@@ -5680,7 +6276,10 @@
 	  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 	  L = CHOLMOD_NAME(analyze) (A, cm);
 	  CHOLMOD_NAME(factorize) (A, L, cm);
-	  rcond = CHOLMOD_NAME(rcond)(L, cm);
+	  if (calc_cond)
+	    rcond = CHOLMOD_NAME(rcond)(L, cm);
+	  else
+	    rcond = 1.;
 	  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 
 	  if (rcond == 0.0)
@@ -5698,7 +6297,10 @@
 		  err = -2;
 
 		  if (sing_handler)
-		    sing_handler (rcond);
+		    {
+		      sing_handler (rcond);
+		      mattype.mark_as_rectangular ();
+		    }
 		  else
 		    (*current_liboctave_error_handler)
 		      ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -5740,7 +6342,8 @@
 	{
 #ifdef HAVE_UMFPACK
 	  Matrix Control, Info;
-	  void *Numeric = factorize (err, rcond, Control, Info, sing_handler);
+	  void *Numeric = factorize (err, rcond, Control, Info,
+				     sing_handler, calc_cond);
 
 	  if (err == 0)
 	    {
@@ -5783,6 +6386,9 @@
 
 	      UMFPACK_ZNAME (free_numeric) (&Numeric);
 	    }
+	  else
+	    mattype.mark_as_rectangular ();
+
 #else
 	  (*current_liboctave_error_handler) ("UMFPACK not installed");
 #endif
@@ -5797,7 +6403,8 @@
 SparseComplexMatrix
 SparseComplexMatrix::fsolve (SparseType &mattype, const SparseComplexMatrix& b,
 			     octave_idx_type& err, double& rcond,
-			     solve_singularity_handler sing_handler) const
+			     solve_singularity_handler sing_handler,
+			     bool calc_cond) const
 {
   SparseComplexMatrix retval;
 
@@ -5916,7 +6523,10 @@
 	  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 	  L = CHOLMOD_NAME(analyze) (A, cm);
 	  CHOLMOD_NAME(factorize) (A, L, cm);
-	  rcond = CHOLMOD_NAME(rcond)(L, cm);
+	  if (calc_cond)
+	    rcond = CHOLMOD_NAME(rcond)(L, cm);
+	  else
+	    rcond = 1.;
 	  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 
 	  if (rcond == 0.0)
@@ -5934,7 +6544,10 @@
 		  err = -2;
 
 		  if (sing_handler)
-		    sing_handler (rcond);
+		    {
+		      sing_handler (rcond);
+		      mattype.mark_as_rectangular ();
+		    }
 		  else
 		    (*current_liboctave_error_handler)
 		      ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
@@ -5982,7 +6595,8 @@
 	{
 #ifdef HAVE_UMFPACK
 	  Matrix Control, Info;
-	  void *Numeric = factorize (err, rcond, Control, Info, sing_handler);
+	  void *Numeric = factorize (err, rcond, Control, Info,
+				     sing_handler, calc_cond);
 
 	  if (err == 0)
 	    {
@@ -5999,7 +6613,7 @@
 
 	      // Take a first guess that the number of non-zero terms
 	      // will be as many as in b
-	      octave_idx_type x_nz = b.nzmax ();
+	      octave_idx_type x_nz = b.nnz ();
 	      octave_idx_type ii = 0;
 	      retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
 
@@ -6072,6 +6686,9 @@
 
 	      UMFPACK_ZNAME (free_numeric) (&Numeric);
 	    }
+	  else
+	    mattype.mark_as_rectangular ();
+
 #else
 	  (*current_liboctave_error_handler) ("UMFPACK not installed");
 #endif
@@ -6111,30 +6728,43 @@
 			    double& rcond, 
 			    solve_singularity_handler sing_handler) const
 {
+  ComplexMatrix retval;
   int typ = mattype.type (false);
 
   if (typ == SparseType::Unknown)
     typ = mattype.type (*this);
 
   if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)
-    return dsolve (mattype, b, err, rcond, sing_handler);
+    retval = dsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper)
-    return utsolve (mattype, b, err, rcond, sing_handler);
+    retval = utsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower)
-    return ltsolve (mattype, b, err, rcond, sing_handler);
+    retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian)
-    return bsolve (mattype, b, err, rcond, sing_handler);
+    retval = bsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Tridiagonal || 
 	   typ == SparseType::Tridiagonal_Hermitian)
-    return trisolve (mattype, b, err, rcond, sing_handler);
+    retval = trisolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Full || typ == SparseType::Hermitian)
-    return fsolve (mattype, b, err, rcond, sing_handler);
-  else
+    retval = fsolve (mattype, b, err, rcond, sing_handler, true);
+  else if (typ != SparseType::Rectangular)
     {
-      (*current_liboctave_error_handler) 
-	("matrix dimension mismatch solution of linear equations");
+      (*current_liboctave_error_handler) ("unknown matrix type");
       return ComplexMatrix ();
     }
+
+  if (mattype.type(false) == SparseType::Rectangular)
+    {
+      rcond = 1.;
+#ifdef USE_QRSOLVE
+      retval = qrsolve (*this, b, err);
+#else
+      retval = dmsolve<ComplexMatrix, SparseComplexMatrix,
+	Matrix> (*this, b, err);
+#endif
+    }
+
+  return retval;
 }
 
 SparseComplexMatrix
@@ -6165,30 +6795,43 @@
 			    octave_idx_type& err, double& rcond,
 			    solve_singularity_handler sing_handler) const
 {
+  SparseComplexMatrix retval;
   int typ = mattype.type (false);
 
   if (typ == SparseType::Unknown)
     typ = mattype.type (*this);
 
   if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)
-    return dsolve (mattype, b, err, rcond, sing_handler);
+    retval = dsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper)
-    return utsolve (mattype, b, err, rcond, sing_handler);
+    retval = utsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower)
-    return ltsolve (mattype, b, err, rcond, sing_handler);
+    retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian)
-    return bsolve (mattype, b, err, rcond, sing_handler);
+    retval = bsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Tridiagonal || 
 	   typ == SparseType::Tridiagonal_Hermitian)
-    return trisolve (mattype, b, err, rcond, sing_handler);
+    retval = trisolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Full || typ == SparseType::Hermitian)
-    return fsolve (mattype, b, err, rcond, sing_handler);
-  else
+    retval = fsolve (mattype, b, err, rcond, sing_handler, true);
+  else if (typ != SparseType::Rectangular)
     {
-      (*current_liboctave_error_handler) 
-	("matrix dimension mismatch solution of linear equations");
+      (*current_liboctave_error_handler) ("unknown matrix type");
       return SparseComplexMatrix ();
     }
+
+  if (mattype.type(false) == SparseType::Rectangular)
+    {
+      rcond = 1.;
+#ifdef USE_QRSOLVE
+      retval = qrsolve (*this, b, err);
+#else
+      retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix,
+	SparseMatrix> (*this, b, err);
+#endif
+    }
+
+  return retval;
 }
 
 ComplexMatrix
@@ -6219,30 +6862,43 @@
 		     octave_idx_type& err, double& rcond, 
 		     solve_singularity_handler sing_handler) const
 {
+  ComplexMatrix retval;
   int typ = mattype.type (false);
 
   if (typ == SparseType::Unknown)
     typ = mattype.type (*this);
 
   if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)
-    return dsolve (mattype, b, err, rcond, sing_handler);
+    retval = dsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper)
-    return utsolve (mattype, b, err, rcond, sing_handler);
+    retval = utsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower)
-    return ltsolve (mattype, b, err, rcond, sing_handler);
+    retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian)
-    return bsolve (mattype, b, err, rcond, sing_handler);
+    retval = bsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Tridiagonal || 
 	   typ == SparseType::Tridiagonal_Hermitian)
-    return trisolve (mattype, b, err, rcond, sing_handler);
+    retval = trisolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Full || typ == SparseType::Hermitian)
-    return fsolve (mattype, b, err, rcond, sing_handler);
-  else
+    retval = fsolve (mattype, b, err, rcond, sing_handler, true);
+  else if (typ != SparseType::Rectangular)
     {
-      (*current_liboctave_error_handler) 
-	("matrix dimension mismatch solution of linear equations");
+      (*current_liboctave_error_handler) ("unknown matrix type");
       return ComplexMatrix ();
     }
+
+  if (mattype.type(false) == SparseType::Rectangular)
+    {
+      rcond = 1.;
+#ifdef USE_QRSOLVE
+      retval = qrsolve (*this, b, err);
+#else
+      retval = dmsolve<ComplexMatrix, SparseComplexMatrix,
+	ComplexMatrix> (*this, b, err);
+#endif
+    }
+
+  return retval;
 }
 
 SparseComplexMatrix
@@ -6274,30 +6930,43 @@
 			    octave_idx_type& err, double& rcond,
 			    solve_singularity_handler sing_handler) const
 {
+  SparseComplexMatrix retval;
   int typ = mattype.type (false);
 
   if (typ == SparseType::Unknown)
     typ = mattype.type (*this);
 
   if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)
-    return dsolve (mattype, b, err, rcond, sing_handler);
+    retval = dsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper)
-    return utsolve (mattype, b, err, rcond, sing_handler);
+    retval = utsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower)
-    return ltsolve (mattype, b, err, rcond, sing_handler);
+    retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian)
-    return bsolve (mattype, b, err, rcond, sing_handler);
+    retval = bsolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Tridiagonal || 
 	   typ == SparseType::Tridiagonal_Hermitian)
-    return trisolve (mattype, b, err, rcond, sing_handler);
+    retval = trisolve (mattype, b, err, rcond, sing_handler, false);
   else if (typ == SparseType::Full || typ == SparseType::Hermitian)
-    return fsolve (mattype, b, err, rcond, sing_handler);
-  else
+    retval = fsolve (mattype, b, err, rcond, sing_handler, true);
+  else if (typ != SparseType::Rectangular)
     {
-      (*current_liboctave_error_handler) 
-	("matrix dimension mismatch solution of linear equations");
+      (*current_liboctave_error_handler) ("unknown matrix type");
       return SparseComplexMatrix ();
     }
+
+  if (mattype.type(false) == SparseType::Rectangular)
+    {
+      rcond = 1.;
+#ifdef USE_QRSOLVE
+      retval = qrsolve (*this, b, err);
+#else
+      retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix,
+	SparseComplexMatrix> (*this, b, err);
+#endif
+    }
+
+  return retval;
 }
 
 ComplexColumnVector
@@ -6543,145 +7212,13 @@
   return solve (tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0));
 }
 
-ComplexMatrix
-SparseComplexMatrix::lssolve (const Matrix& b) const
-{
-  octave_idx_type info;
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-ComplexMatrix
-SparseComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info) const
-{
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-ComplexMatrix
-SparseComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type&) const
-{
-  return qrsolve (*this, b, info);
-}
-
-SparseComplexMatrix
-SparseComplexMatrix::lssolve (const SparseMatrix& b) const
-{
-  octave_idx_type info;
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-SparseComplexMatrix
-SparseComplexMatrix::lssolve (const SparseMatrix& b, octave_idx_type& info) const
-{
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-SparseComplexMatrix
-SparseComplexMatrix::lssolve (const SparseMatrix& b, octave_idx_type& info, 
-			      octave_idx_type&) const
-{
-  return qrsolve (*this, b, info);
-}
-
-ComplexMatrix
-SparseComplexMatrix::lssolve (const ComplexMatrix& b) const
-{
-  octave_idx_type info;
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-ComplexMatrix
-SparseComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info) const
-{
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-ComplexMatrix
-SparseComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, 
-			      octave_idx_type&) const
-{
-  return qrsolve (*this, b, info);
-}
-
-SparseComplexMatrix
-SparseComplexMatrix::lssolve (const SparseComplexMatrix& b) const
-{
-  octave_idx_type info;
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-SparseComplexMatrix
-SparseComplexMatrix::lssolve (const SparseComplexMatrix& b, octave_idx_type& info) const
-{
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-SparseComplexMatrix
-SparseComplexMatrix::lssolve (const SparseComplexMatrix& b, octave_idx_type& info, 
-			      octave_idx_type&) const
-{
-  return qrsolve (*this, b, info);
-}
-
-ComplexColumnVector
-SparseComplexMatrix::lssolve (const ColumnVector& b) const
-{
-  octave_idx_type info;
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-ComplexColumnVector
-SparseComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info) const
-{
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-ComplexColumnVector
-SparseComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const
-{
-  Matrix tmp (b);
-  return lssolve (tmp, info, rank).column (static_cast<octave_idx_type> (0));
-}
-
-ComplexColumnVector
-SparseComplexMatrix::lssolve (const ComplexColumnVector& b) const
-{
-  octave_idx_type info;
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-ComplexColumnVector
-SparseComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info) const
-{
-  octave_idx_type rank;
-  return lssolve (b, info, rank);
-}
-
-ComplexColumnVector
-SparseComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info,
-			octave_idx_type& rank) const
-{
-  ComplexMatrix tmp (b);
-  return lssolve (tmp, info, rank).column (static_cast<octave_idx_type> (0));
-}
-
 // unary operations
 SparseBoolMatrix
 SparseComplexMatrix::operator ! (void) const
 {
   octave_idx_type nr = rows ();
   octave_idx_type nc = cols ();
-  octave_idx_type nz1 = nzmax ();
+  octave_idx_type nz1 = nnz ();
   octave_idx_type nz2 = nr*nc - nz1;
    
   SparseBoolMatrix r (nr, nc, nz2);
@@ -6755,7 +7292,7 @@
 {
   octave_idx_type nr = rows ();
   octave_idx_type nc = cols ();
-  octave_idx_type nz = nzmax ();
+  octave_idx_type nz = nnz ();
   bool f_zero = (f(0.0) == 0.0);
 
   // Count number of non-zero elements
@@ -6805,7 +7342,7 @@
 {
   octave_idx_type nr = rows ();
   octave_idx_type nc = cols ();
-  octave_idx_type nz = nzmax ();
+  octave_idx_type nz = nnz ();
   bool f_zero = (f(0.0) == 0.0);
 
   // Count number of non-zero elements
@@ -6855,7 +7392,7 @@
 {
   octave_idx_type nr = rows ();
   octave_idx_type nc = cols ();
-  octave_idx_type nz = nzmax ();
+  octave_idx_type nz = nnz ();
   bool f_zero = f(0.0);
 
   // Count number of non-zero elements
@@ -6910,7 +7447,7 @@
 bool
 SparseComplexMatrix::any_element_is_inf_or_nan (void) const
 {
-  octave_idx_type nel = nzmax ();
+  octave_idx_type nel = nnz ();
 
   for (octave_idx_type i = 0; i < nel; i++)
     {
@@ -6927,7 +7464,7 @@
 bool
 SparseComplexMatrix::all_elements_are_real (void) const
 {
-  octave_idx_type nel = nzmax ();
+  octave_idx_type nel = nnz ();
 
   for (octave_idx_type i = 0; i < nel; i++)
     {
@@ -6947,7 +7484,7 @@
 bool
 SparseComplexMatrix::all_integers (double& max_val, double& min_val) const
 {
-  octave_idx_type nel = nzmax ();
+  octave_idx_type nel = nnz ();
 
   if (nel == 0)
     return false;
@@ -6984,7 +7521,7 @@
 bool
 SparseComplexMatrix::too_large_for_float (void) const
 {
-  octave_idx_type nel = nzmax ();
+  octave_idx_type nel = nnz ();
 
   for (octave_idx_type i = 0; i < nel; i++)
     {
@@ -7062,7 +7599,7 @@
 
 SparseMatrix SparseComplexMatrix::abs (void) const
 {
-  octave_idx_type nz = nzmax ();
+  octave_idx_type nz = nnz ();
   octave_idx_type nc = cols ();
 
   SparseMatrix retval (rows(), nc, nz);
@@ -7237,74 +7774,55 @@
 SparseComplexMatrix
 operator * (const SparseComplexMatrix& m, const SparseMatrix& a)
 {
-  SparseComplexMatrix tmp (a);
-  return m * tmp;
+  SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, double);
 }
 
 SparseComplexMatrix
 operator * (const SparseMatrix& m, const SparseComplexMatrix& a)
 {
-  SparseComplexMatrix tmp (m);
-  return tmp * a;
+  SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, Complex);
 }
 
 SparseComplexMatrix
 operator * (const SparseComplexMatrix& m, const SparseComplexMatrix& a)
 {
-#ifdef HAVE_SPARSE_BLAS
-  // XXX FIXME XXX Isn't there a sparse BLAS ??
-#else
-  // Use Andy's sparse matrix multiply function
-  SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex);
-#endif
+  SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, Complex);
 }
 
 ComplexMatrix
 operator * (const ComplexMatrix& m, const SparseMatrix& a)
 {
-  SparseComplexMatrix tmp (a);
-  return m * tmp;
+  FULL_SPARSE_MUL (ComplexMatrix, double, Complex (0.,0.));
 }
 
 ComplexMatrix
 operator * (const Matrix& m, const SparseComplexMatrix& a)
 {
-  ComplexMatrix tmp (m);
-  return tmp * a;
+  FULL_SPARSE_MUL (ComplexMatrix, Complex, Complex (0.,0.));
 }
 
 ComplexMatrix
 operator * (const ComplexMatrix& m, const SparseComplexMatrix& a)
 {
-#ifdef HAVE_SPARSE_BLAS
-  // XXX FIXME XXX Isn't there a sparse BLAS ??
-#else
-  FULL_SPARSE_MUL (ComplexMatrix, Complex);
-#endif
+  FULL_SPARSE_MUL (ComplexMatrix, Complex, Complex (0.,0.));
 }
 
 ComplexMatrix
 operator * (const SparseComplexMatrix& m, const Matrix& a)
 {
-  ComplexMatrix tmp (a);
-  return m * tmp;
+  SPARSE_FULL_MUL (ComplexMatrix, double, Complex (0.,0.));
 }
 
 ComplexMatrix
 operator * (const SparseMatrix& m, const ComplexMatrix& a)
 {
-  SparseComplexMatrix tmp (m);
-  return tmp * a;
+  SPARSE_FULL_MUL (ComplexMatrix, Complex, Complex (0.,0.));
 }
 
 ComplexMatrix
 operator * (const SparseComplexMatrix& m, const ComplexMatrix& a)
 {
-#ifdef HAVE_SPARSE_BLAS
-  // XXX FIXME XXX Isn't there a sparse BLAS ??
-#else
-  SPARSE_FULL_MUL (ComplexMatrix, Complex);
-#endif
+  SPARSE_FULL_MUL (ComplexMatrix, Complex, Complex (0.,0.));
 }
 
 // XXX FIXME XXX -- it would be nice to share code among the min/max
@@ -7357,14 +7875,14 @@
       octave_idx_type b_nr = b.rows ();
       octave_idx_type b_nc = b.cols ();
 
-      if (a_nr == 0 || b_nc == 0 || a.nzmax () == 0 || b.nzmax () == 0)
+      if (a_nr == 0 || b_nc == 0 || a.nnz () == 0 || b.nnz () == 0)
 	return SparseComplexMatrix (a_nr, a_nc);
 
       if (a_nr != b_nr || a_nc != b_nc)
 	gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
       else
 	{
-	  r = SparseComplexMatrix (a_nr, a_nc, (a.nzmax () + b.nzmax ()));
+	  r = SparseComplexMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
        
 	  octave_idx_type jx = 0;
 	  r.cidx (0) = 0;
@@ -7479,16 +7997,16 @@
 
       if (a_nr == 0 || b_nc == 0)
 	return SparseComplexMatrix (a_nr, a_nc);
-      if (a.nzmax () == 0)
+      if (a.nnz () == 0)
 	return SparseComplexMatrix (b);
-      if (b.nzmax () == 0)
+      if (b.nnz () == 0)
 	return SparseComplexMatrix (a);
 
       if (a_nr != b_nr || a_nc != b_nc)
 	gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
       else
 	{
-	  r = SparseComplexMatrix (a_nr, a_nc, (a.nzmax () + b.nzmax ()));
+	  r = SparseComplexMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
        
 	  octave_idx_type jx = 0;
 	  r.cidx (0) = 0;