octave-nkf: liboctave/CSparse.cc comparison

comparison liboctave/CSparse.cc @ 5681:233d98d95659

[project @ 2006-03-16 17:48:55 by dbateman]

author	dbateman
date	Thu, 16 Mar 2006 17:48:56 +0000
parents	69a4f320d95a
children	2fe20065a545

comparison

equal deleted inserted replaced

-:cc6a965ae4ca
+:233d98d95659
 #include "SparseCmplxCHOL.h"
 #include "SparseCmplxQR.h"
 #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"
 {
 F77_RET_T
 F77_FUNC (zgbtrf, ZGBTRF) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
 			     F77_CHAR_ARG_LEN_DECL);
 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
 F77_FUNC (zgttrf, ZGTTRF) (const octave_idx_type&, Complex*, Complex*, Complex*,
 			     Complex*, octave_idx_type*, octave_idx_type&);
 F77_FUNC (zgtsv, ZGTSV) (const octave_idx_type&, const octave_idx_type&, Complex*, Complex*,
 			   Complex*, Complex*, const octave_idx_type&, octave_idx_type&);
 }
 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);
 }
 }
 bool
 SparseComplexMatrix::operator == (const SparseComplexMatrix& a) const
 {
 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;
 for (octave_idx_type i = 0; i < nc + 1; i++)
 SparseComplexMatrix
 SparseComplexMatrix::hermitian (void) const
 {
 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);
 for (octave_idx_type i = 0; i < nr; i++)
 w[i] = 0;
 SparseComplexMatrix
 conj (const SparseComplexMatrix& a)
 {
 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++)
 retval.cidx (i) = a.cidx (i);
 		}
 	    }
 	  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);
 	      for (octave_idx_type i = 0; i < nr; i++)
 	      retval.xcidx(nr) = cx;
 	      retval.maybe_compress ();
 	    }
 	  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);
 	      OCTAVE_LOCAL_BUFFER (Complex, work, nr);
 return retval;
 }
 ComplexMatrix
 SparseComplexMatrix::dsolve (SparseType &mattype, const Matrix& b,
-			     octave_idx_type& err,
+			     octave_idx_type& err, double& rcond,
-			     double& rcond, solve_singularity_handler) const
+			     solve_singularity_handler, bool calc_cond) const
 {
 ComplexMatrix retval;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	    for (octave_idx_type j = 0; j < b.cols(); j++)
 	      for (octave_idx_type k = 0; k < nc; k++)
 		for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
 		  retval(k,j) = b(ridx(i),j) / data (i);
-	  double dmax = 0., dmin = octave_Inf;
+	  if (calc_cond)
-	  for (octave_idx_type i = 0; i < nm; i++)
+	    {
-	    {
+	      double dmax = 0., dmin = octave_Inf;
-	      double tmp = std::abs(data(i));
+	      for (octave_idx_type i = 0; i < nm; i++)
-	      if (tmp > dmax)
+		{
-		dmax = tmp;
+		  double tmp = std::abs(data(i));
-	      if (tmp < dmin)
+		  if (tmp > dmax)
-		dmin = tmp;
+		    dmax = tmp;
-	    }
+		  if (tmp < dmin)
-	  rcond = dmin / dmax;
+		    dmin = tmp;
+		}
+	      rcond = dmin / dmax;
+	    }
+	  else
+	    rcond = 1.0;
 	}
 else
 	(*current_liboctave_error_handler) ("incorrect matrix type");
 }
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (typ == SparseType::Diagonal ||
 	  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;
 	  octave_idx_type ii = 0;
 	  if (typ == SparseType::Diagonal)
 	    for (octave_idx_type j = 0; j < b.cols(); j++)
 	      {
 		for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
 		  {
+		    if (b.ridx(i) >= nm)
+		      break;
 		    retval.xridx (ii) = b.ridx(i);
 		    retval.xdata (ii++) = b.data(i) / data (b.ridx (i));
 		  }
 		retval.xcidx(j+1) = ii;
 	      }
 			}
 		    }
 		retval.xcidx(j+1) = ii;
 	      }
-	  double dmax = 0., dmin = octave_Inf;
+	  if (calc_cond)
-	  for (octave_idx_type i = 0; i < nm; i++)
+	    {
-	    {
+	      double dmax = 0., dmin = octave_Inf;
-	      double tmp = std::abs(data(i));
+	      for (octave_idx_type i = 0; i < nm; i++)
-	      if (tmp > dmax)
+		{
-		dmax = tmp;
+		  double tmp = std::abs(data(i));
-	      if (tmp < dmin)
+		  if (tmp > dmax)
-		dmin = tmp;
+		    dmax = tmp;
-	    }
+		  if (tmp < dmin)
-	  rcond = dmin / dmax;
+		    dmin = tmp;
+		}
+	      rcond = dmin / dmax;
+	    }
+	  else
+	    rcond = 1.0;
 	}
 else
 	(*current_liboctave_error_handler) ("incorrect matrix type");
 }
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	    for (octave_idx_type j = 0; j < b.cols(); j++)
 	      for (octave_idx_type k = 0; k < nc; k++)
 		for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
 		  retval(k,j) = b(ridx(i),j) / data (i);
-	  double dmax = 0., dmin = octave_Inf;
+	  if (calc_cond)
-	  for (octave_idx_type i = 0; i < nr; i++)
+	    {
-	    {
+	      double dmax = 0., dmin = octave_Inf;
-	      double tmp = std::abs(data(i));
+	      for (octave_idx_type i = 0; i < nr; i++)
-	      if (tmp > dmax)
+		{
-		dmax = tmp;
+		  double tmp = std::abs(data(i));
-	      if (tmp < dmin)
+		  if (tmp > dmax)
-		dmin = tmp;
+		    dmax = tmp;
-	    }
+		  if (tmp < dmin)
-	  rcond = dmin / dmax;
+		    dmin = tmp;
+		}
+	      rcond = dmin / dmax;
+	    }
+	  else
+	    rcond = 1.0;
 	}
 else
 	(*current_liboctave_error_handler) ("incorrect matrix type");
 }
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (typ == SparseType::Diagonal ||
 	  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;
 	  octave_idx_type ii = 0;
 	  if (typ == SparseType::Diagonal)
 	    for (octave_idx_type j = 0; j < b.cols(); j++)
 	      {
 		for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
 		  {
+		    if (b.ridx(i) >= nm)
+		      break;
 		    retval.xridx (ii) = b.ridx(i);
 		    retval.xdata (ii++) = b.data(i) / data (b.ridx (i));
 		  }
 		retval.xcidx(j+1) = ii;
 	      }
 			}
 		    }
 		retval.xcidx(j+1) = ii;
 	      }
-	  double dmax = 0., dmin = octave_Inf;
+	  if (calc_cond)
-	  for (octave_idx_type i = 0; i < nm; i++)
+	    {
-	    {
+	      double dmax = 0., dmin = octave_Inf;
-	      double tmp = std::abs(data(i));
+	      for (octave_idx_type i = 0; i < nm; i++)
-	      if (tmp > dmax)
+		{
-		dmax = tmp;
+		  double tmp = std::abs(data(i));
-	      if (tmp < dmin)
+		  if (tmp > dmax)
-		dmin = tmp;
+		    dmax = tmp;
-	    }
+		  if (tmp < dmin)
-	  rcond = dmin / dmax;
+		    dmin = tmp;
+		}
+	      rcond = dmin / dmax;
+	    }
+	  else
+	    rcond = 1.0;
 	}
 else
 	(*current_liboctave_error_handler) ("incorrect matrix type");
 }
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	  typ == SparseType::Upper)
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
 	  octave_idx_type b_nc = b.cols ();
-	  rcond = 0.;
+	  rcond = 1.;
-	  // Calculate the 1-norm of matrix for rcond calculation
+	  if (calc_cond)
-	  for (octave_idx_type j = 0; j < nc; j++)
+	    {
-	    {
+	      // Calculate the 1-norm of matrix for rcond calculation
-	      double atmp = 0.;
+	      for (octave_idx_type j = 0; j < nc; j++)
-	      for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		{
-		atmp += std::abs(data(i));
+		  double atmp = 0.;
-	      if (atmp > anorm)
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
-		anorm = atmp;
+		    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++)
 		{
 		  for (octave_idx_type i = 0; i < nr; i++)
 		    work[i] = b(i,j);
 		    {
 		      octave_idx_type kidx = perm[k];
 		      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;
 			    }
 		  for (octave_idx_type i = 0; i < nc; i++)
 		    retval (perm[i], j) = work[i];
 		}
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k >= 0; k--)
+		      work[j] = 1.;
-		    {
-		      octave_idx_type iidx = perm[k];
+		      for (octave_idx_type k = j; k >= 0; k--)
-		      if (work[k] != 0.)
 			{
-			  Complex tmp = work[k] / data(cidx(iidx+1)-1);
+			  octave_idx_type iidx = perm[k];
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(iidx);
+			  if (work[k] != 0.)
-			       i < cidx(iidx+1)-1; i++)
 			    {
-			      octave_idx_type idx2 = ridx(i);
+			      Complex tmp = work[k] / data(cidx(iidx+1)-1);
-			      work[idx2] = work[idx2] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
-		  for (octave_idx_type i = 0; i < j+1; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    ainvnorm = atmp;
+		    }
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
 	    {
 	      OCTAVE_LOCAL_BUFFER (Complex, work, nm);
 		  for (octave_idx_type k = nc-1; k >= 0; k--)
 		    {
 		      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;
 			    }
 		  for (octave_idx_type i = 0; i < nc; i++)
 		    retval.xelem (i, j) = work[i];
 		}
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k >= 0; k--)
+		      work[j] = 1.;
-		    {
-		      if (work[k] != 0.)
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex tmp = work[k] / data(cidx(k+1)-1);
+			  if (work[k] != 0.)
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
 			    {
-			      octave_idx_type iidx = ridx(i);
+			      Complex tmp = work[k] / data(cidx(k+1)-1);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
-		  for (octave_idx_type i = 0; i < j+1; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    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",
 		   rcond);
 	    }
 	  if (rcond_plus_one == 1.0 || xisnan (rcond))
 	    {
 	      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",
 		   rcond);
 	    }
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (typ == SparseType::Permuted_Upper ||
 	  typ == SparseType::Upper)
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
-	  rcond = 0.;
+	  rcond = 1.;
-	  // Calculate the 1-norm of matrix for rcond calculation
+	  if (calc_cond)
-	  for (octave_idx_type j = 0; j < nc; j++)
+	    {
-	    {
+	      // Calculate the 1-norm of matrix for rcond calculation
-	      double atmp = 0.;
+	      for (octave_idx_type j = 0; j < nc; j++)
-	      for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		{
-		atmp += std::abs(data(i));
+		  double atmp = 0.;
-	      if (atmp > anorm)
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
-		anorm = atmp;
+		    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;
 	  octave_idx_type x_nz = b_nz;
 		    {
 		      octave_idx_type kidx = perm[k];
 		      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;
 			    }
 		  retval.xcidx(j+1) = ii;
 		}
 	      retval.maybe_compress ();
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k >= 0; k--)
+		      work[j] = 1.;
-		    {
-		      octave_idx_type iidx = perm[k];
+		      for (octave_idx_type k = j; k >= 0; k--)
-		      if (work[k] != 0.)
 			{
-			  Complex tmp = work[k] / data(cidx(iidx+1)-1);
+			  octave_idx_type iidx = perm[k];
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(iidx);
+			  if (work[k] != 0.)
-			       i < cidx(iidx+1)-1; i++)
 			    {
-			      octave_idx_type idx2 = ridx(i);
+			      Complex tmp = work[k] / data(cidx(iidx+1)-1);
-			      work[idx2] = work[idx2] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
-		  for (octave_idx_type i = 0; i < j+1; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    ainvnorm = atmp;
+		    }
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
 	    {
 	      OCTAVE_LOCAL_BUFFER (Complex, work, nm);
 		  for (octave_idx_type k = nc-1; k >= 0; k--)
 		    {
 		      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;
 			    }
 		  retval.xcidx(j+1) = ii;
 		}
 	      retval.maybe_compress ();
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k >= 0; k--)
+		      work[j] = 1.;
-		    {
-		      if (work[k] != 0.)
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex tmp = work[k] / data(cidx(k+1)-1);
+			  if (work[k] != 0.)
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
 			    {
-			      octave_idx_type iidx = ridx(i);
+			      Complex tmp = work[k] / data(cidx(k+1)-1);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
-		  for (octave_idx_type i = 0; i < j+1; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    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",
 		   rcond);
 	    }
 	  if (rcond_plus_one == 1.0 || xisnan (rcond))
 	    {
 	      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",
 		   rcond);
 	    }
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	  typ == SparseType::Upper)
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
 	  octave_idx_type b_nc = b.cols ();
-	  rcond = 0.;
+	  rcond = 1.;
-	  // Calculate the 1-norm of matrix for rcond calculation
+	  if (calc_cond)
-	  for (octave_idx_type j = 0; j < nc; j++)
+	    {
-	    {
+	      // Calculate the 1-norm of matrix for rcond calculation
-	      double atmp = 0.;
+	      for (octave_idx_type j = 0; j < nc; j++)
-	      for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		{
-		atmp += std::abs(data(i));
+		  double atmp = 0.;
-	      if (atmp > anorm)
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
-		anorm = atmp;
+		    atmp += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
 	    }
 	  if (typ == SparseType::Permuted_Upper)
 	    {
 	      retval.resize (nc, b_nc);
 		    {
 		      octave_idx_type kidx = perm[k];
 		      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;
 			    }
 		  for (octave_idx_type i = 0; i < nc; i++)
 		    retval (perm[i], j) = work[i];
 		}
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k >= 0; k--)
+		      work[j] = 1.;
-		    {
-		      octave_idx_type iidx = perm[k];
+		      for (octave_idx_type k = j; k >= 0; k--)
-		      if (work[k] != 0.)
 			{
-			  Complex tmp = work[k] / data(cidx(iidx+1)-1);
+			  octave_idx_type iidx = perm[k];
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(iidx);
+			  if (work[k] != 0.)
-			       i < cidx(iidx+1)-1; i++)
 			    {
-			      octave_idx_type idx2 = ridx(i);
+			      Complex tmp = work[k] / data(cidx(iidx+1)-1);
-			      work[idx2] = work[idx2] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
-		  for (octave_idx_type i = 0; i < j+1; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    ainvnorm = atmp;
+		    }
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
 	    {
 	      OCTAVE_LOCAL_BUFFER (Complex, work, nm);
 		  for (octave_idx_type k = nc-1; k >= 0; k--)
 		    {
 		      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;
 			    }
 		  for (octave_idx_type i = 0; i < nc; i++)
 		    retval.xelem (i, j) = work[i];
 		}
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k >= 0; k--)
+		      work[j] = 1.;
-		    {
-		      if (work[k] != 0.)
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex tmp = work[k] / data(cidx(k+1)-1);
+			  if (work[k] != 0.)
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
 			    {
-			      octave_idx_type iidx = ridx(i);
+			      Complex tmp = work[k] / data(cidx(k+1)-1);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
-		  for (octave_idx_type i = 0; i < j+1; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    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",
 		   rcond);
 	    }
 	  if (rcond_plus_one == 1.0 || xisnan (rcond))
 	    {
 	      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",
 		   rcond);
 	    }
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (typ == SparseType::Permuted_Upper ||
 	  typ == SparseType::Upper)
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
-	  rcond = 0.;
+	  rcond = 1.;
-	  // Calculate the 1-norm of matrix for rcond calculation
+	  if (calc_cond)
-	  for (octave_idx_type j = 0; j < nc; j++)
+	    {
-	    {
+	      // Calculate the 1-norm of matrix for rcond calculation
-	      double atmp = 0.;
+	      for (octave_idx_type j = 0; j < nc; j++)
-	      for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		{
-		atmp += std::abs(data(i));
+		  double atmp = 0.;
-	      if (atmp > anorm)
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
-		anorm = atmp;
+		    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;
 	  octave_idx_type x_nz = b_nz;
 		    {
 		      octave_idx_type kidx = perm[k];
 		      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;
 			    }
 		  retval.xcidx(j+1) = ii;
 		}
 	      retval.maybe_compress ();
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k >= 0; k--)
+		      work[j] = 1.;
-		    {
-		      octave_idx_type iidx = perm[k];
+		      for (octave_idx_type k = j; k >= 0; k--)
-		      if (work[k] != 0.)
 			{
-			  Complex tmp = work[k] / data(cidx(iidx+1)-1);
+			  octave_idx_type iidx = perm[k];
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(iidx);
+			  if (work[k] != 0.)
-			       i < cidx(iidx+1)-1; i++)
 			    {
-			      octave_idx_type idx2 = ridx(i);
+			      Complex tmp = work[k] / data(cidx(iidx+1)-1);
-			      work[idx2] = work[idx2] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
-		  for (octave_idx_type i = 0; i < j+1; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    ainvnorm = atmp;
+		    }
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
 	    {
 	      OCTAVE_LOCAL_BUFFER (Complex, work, nm);
 		  for (octave_idx_type k = nr-1; k >= 0; k--)
 		    {
 		      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;
 			    }
 		  retval.xcidx(j+1) = ii;
 		}
 	      retval.maybe_compress ();
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k >= 0; k--)
+		      work[j] = 1.;
-		    {
-		      if (work[k] != 0.)
+		      for (octave_idx_type k = j; k >= 0; k--)
 			{
-			  Complex tmp = work[k] / data(cidx(k+1)-1);
+			  if (work[k] != 0.)
-			  work[k] = tmp;
-			  for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
 			    {
-			      octave_idx_type iidx = ridx(i);
+			      Complex tmp = work[k] / data(cidx(k+1)-1);
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = 0; i < j+1; i++)
-		  for (octave_idx_type i = 0; i < j+1; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    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",
 		   rcond);
 	    }
 	  if (rcond_plus_one == 1.0 || xisnan (rcond))
 	    {
 	      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",
 		   rcond);
 	    }
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	  typ == SparseType::Lower)
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
 	  octave_idx_type b_nc = b.cols ();
-	  rcond = 0.;
+	  rcond = 1.;
-	  // Calculate the 1-norm of matrix for rcond calculation
+	  if (calc_cond)
-	  for (octave_idx_type j = 0; j < nc; j++)
+	    {
-	    {
+	      // Calculate the 1-norm of matrix for rcond calculation
-	      double atmp = 0.;
+	      for (octave_idx_type j = 0; j < nc; j++)
-	      for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		{
-		atmp += std::abs(data(i));
+		  double atmp = 0.;
-	      if (atmp > anorm)
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
-		anorm = atmp;
+		    atmp += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
 	    }
 	  if (typ == SparseType::Permuted_Lower)
 	    {
 	      retval.resize (nc, b_nc);
 			      {
 				minr = perm[ridx(i)];
 				mini = i;
 			      }
-			  if (minr != k)
+			  if (minr != k || data (mini) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
 			    }
 		  for (octave_idx_type i = 0; i < nc; i++)
 		    retval (i, j) = work[i];
 		}
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = 0; k < nc; k++)
+		      work[j] = 1.;
-		    {
-		      if (work[k] != 0.)
+		      for (octave_idx_type k = 0; k < nc; k++)
 			{
-			  octave_idx_type minr = nr;
+			  if (work[k] != 0.)
-			  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)
+			      octave_idx_type minr = nr;
-				continue;
+			      octave_idx_type mini = 0;
-			      octave_idx_type iidx = perm[ridx(i)];
+			      for (octave_idx_type i = cidx(k);
-			      work[iidx] = work[iidx] - tmp * data(i);
+				   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;
-		  double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
-		  for (octave_idx_type i = j; i < nc; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    ainvnorm = atmp;
+		    }
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
 	    {
 	      OCTAVE_LOCAL_BUFFER (Complex, work, nm);
 		    work[i] = 0.;
 		  for (octave_idx_type k = 0; k < nc; k++)
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k)) != k)
+			  if (ridx(cidx(k)) != k ||
+			      data(cidx(k)) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
 			    }
 		    }
 		  for (octave_idx_type i = 0; i < nc; i++)
 		    retval.xelem (i, j) = work[i];
 		}
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k < nc; k++)
+		      work[j] = 1.;
-		    {
+		      for (octave_idx_type k = j; k < nc; k++)
-		      if (work[k] != 0.)
 			{
-			  Complex tmp = work[k] / data(cidx(k));
-			  work[k] = tmp;
+			  if (work[k] != 0.)
-			  for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
 			    {
-			      octave_idx_type iidx = ridx(i);
+			      Complex tmp = work[k] / data(cidx(k));
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
-		  for (octave_idx_type i = j; i < nc; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    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",
 		   rcond);
 	    }
 	  if (rcond_plus_one == 1.0 || xisnan (rcond))
 	    {
 	      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",
 		   rcond);
 	    }
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (typ == SparseType::Permuted_Lower ||
 	  typ == SparseType::Lower)
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
-	  rcond = 0.;
+	  rcond = 1.;
-	  // Calculate the 1-norm of matrix for rcond calculation
+	  if (calc_cond)
-	  for (octave_idx_type j = 0; j < nc; j++)
+	    {
-	    {
+	      // Calculate the 1-norm of matrix for rcond calculation
-	      double atmp = 0.;
+	      for (octave_idx_type j = 0; j < nc; j++)
-	      for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		{
-		atmp += std::abs(data(i));
+		  double atmp = 0.;
-	      if (atmp > anorm)
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
-		anorm = atmp;
+		    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;
 	  octave_idx_type x_nz = b_nz;
 			      {
 				minr = perm[ridx(i)];
 				mini = i;
 			      }
-			  if (minr != k)
+			  if (minr != k || data (mini) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
 			    }
 		  retval.xcidx(j+1) = ii;
 		}
 	      retval.maybe_compress ();
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = 0; k < nc; k++)
+		      work[j] = 1.;
-		    {
-		      if (work[k] != 0.)
+		      for (octave_idx_type k = 0; k < nc; k++)
 			{
-			  octave_idx_type minr = nr;
+			  if (work[k] != 0.)
-			  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)
+			      octave_idx_type minr = nr;
-				continue;
+			      octave_idx_type mini = 0;
-			      octave_idx_type iidx = perm[ridx(i)];
+			      for (octave_idx_type i = cidx(k);
-			      work[iidx] = work[iidx] - tmp * data(i);
+				   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;
-		  double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
-		  for (octave_idx_type i = j; i < nc; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    ainvnorm = atmp;
+		    }
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
 	    {
 	      OCTAVE_LOCAL_BUFFER (Complex, work, nm);
 		  for (octave_idx_type k = 0; k < nc; k++)
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k)) != k)
+			  if (ridx(cidx(k)) != k ||
+			      data(cidx(k)) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
 			    }
 		  retval.xcidx(j+1) = ii;
 		}
 	      retval.maybe_compress ();
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k < nc; k++)
+		      work[j] = 1.;
-		    {
+		      for (octave_idx_type k = j; k < nc; k++)
-		      if (work[k] != 0.)
 			{
-			  Complex tmp = work[k] / data(cidx(k));
-			  work[k] = tmp;
+			  if (work[k] != 0.)
-			  for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
 			    {
-			      octave_idx_type iidx = ridx(i);
+			      Complex tmp = work[k] / data(cidx(k));
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
-		  for (octave_idx_type i = j; i < nc; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    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",
 		   rcond);
 	    }
 	  if (rcond_plus_one == 1.0 || xisnan (rcond))
 	    {
 	      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",
 		   rcond);
 	    }
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	  typ == SparseType::Lower)
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
 	  octave_idx_type b_nc = b.cols ();
-	  rcond = 0.;
+	  rcond = 1.;
-	  // Calculate the 1-norm of matrix for rcond calculation
+	  if (calc_cond)
-	  for (octave_idx_type j = 0; j < nc; j++)
+	    {
-	    {
+	      // Calculate the 1-norm of matrix for rcond calculation
-	      double atmp = 0.;
+	      for (octave_idx_type j = 0; j < nc; j++)
-	      for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		{
-		atmp += std::abs(data(i));
+		  double atmp = 0.;
-	      if (atmp > anorm)
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
-		anorm = atmp;
+		    atmp += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
 	    }
 	  if (typ == SparseType::Permuted_Lower)
 	    {
 	      retval.resize (nc, b_nc);
 			      {
 				minr = perm[ridx(i)];
 				mini = i;
 			      }
-			  if (minr != k)
+			  if (minr != k || data (mini) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
 			    }
 		  for (octave_idx_type i = 0; i < nc; i++)
 		    retval (i, j) = work[i];
 		}
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = 0; k < nc; k++)
+		      work[j] = 1.;
-		    {
-		      if (work[k] != 0.)
+		      for (octave_idx_type k = 0; k < nc; k++)
 			{
-			  octave_idx_type minr = nr;
+			  if (work[k] != 0.)
-			  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)
+			      octave_idx_type minr = nr;
-				continue;
+			      octave_idx_type mini = 0;
-			      octave_idx_type iidx = perm[ridx(i)];
+			      for (octave_idx_type i = cidx(k);
-			      work[iidx] = work[iidx] - tmp * data(i);
+				   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;
-		  double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
-		  for (octave_idx_type i = j; i < nc; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    ainvnorm = atmp;
+		    }
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
 	    {
 	      OCTAVE_LOCAL_BUFFER (Complex, work, nm);
 		  for (octave_idx_type k = 0; k < nc; k++)
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k)) != k)
+			  if (ridx(cidx(k)) != k ||
+			      data(cidx(k)) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
 			    }
 		  for (octave_idx_type i = 0; i < nc; i++)
 		    retval.xelem (i, j) = work[i];
 		}
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k < nc; k++)
+		      work[j] = 1.;
-		    {
+		      for (octave_idx_type k = j; k < nc; k++)
-		      if (work[k] != 0.)
 			{
-			  Complex tmp = work[k] / data(cidx(k));
-			  work[k] = tmp;
+			  if (work[k] != 0.)
-			  for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
 			    {
-			      octave_idx_type iidx = ridx(i);
+			      Complex tmp = work[k] / data(cidx(k));
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
-		  for (octave_idx_type i = j; i < nc; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    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",
 		   rcond);
 	    }
 	  if (rcond_plus_one == 1.0 || xisnan (rcond))
 	    {
 	      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",
 		   rcond);
 	    }
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 if (typ == SparseType::Permuted_Lower ||
 	  typ == SparseType::Lower)
 	{
 	  double anorm = 0.;
 	  double ainvnorm = 0.;
-	  rcond = 0.;
+	  rcond = 1.;
-	  // Calculate the 1-norm of matrix for rcond calculation
+	  if (calc_cond)
-	  for (octave_idx_type j = 0; j < nc; j++)
+	    {
-	    {
+	      // Calculate the 1-norm of matrix for rcond calculation
-	      double atmp = 0.;
+	      for (octave_idx_type j = 0; j < nc; j++)
-	      for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		{
-		atmp += std::abs(data(i));
+		  double atmp = 0.;
-	      if (atmp > anorm)
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
-		anorm = atmp;
+		    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;
 	  octave_idx_type x_nz = b_nz;
 			      {
 				minr = perm[ridx(i)];
 				mini = i;
 			      }
-			  if (minr != k)
+			  if (minr != k || data (mini) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
 			    }
 		  retval.xcidx(j+1) = ii;
 		}
 	      retval.maybe_compress ();
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = 0; k < nc; k++)
+		      work[j] = 1.;
-		    {
-		      if (work[k] != 0.)
+		      for (octave_idx_type k = 0; k < nc; k++)
 			{
-			  octave_idx_type minr = nr;
+			  if (work[k] != 0.)
-			  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)
+			      octave_idx_type minr = nr;
-				continue;
+			      octave_idx_type mini = 0;
-			      octave_idx_type iidx = perm[ridx(i)];
+			      for (octave_idx_type i = cidx(k);
-			      work[iidx] = work[iidx] - tmp * data(i);
+				   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;
-		  double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
-		  for (octave_idx_type i = j; i < nc; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    ainvnorm = atmp;
+		    }
+		  rcond = 1. / ainvnorm / anorm;
 		}
 	    }
 	  else
 	    {
 	      OCTAVE_LOCAL_BUFFER (Complex, work, nm);
 		  for (octave_idx_type k = 0; k < nc; k++)
 		    {
 		      if (work[k] != 0.)
 			{
-			  if (ridx(cidx(k)) != k)
+			  if (ridx(cidx(k)) != k ||
+			      data(cidx(k)) == 0.)
 			    {
 			      err = -2;
 			      goto triangular_error;
 			    }
 		  retval.xcidx(j+1) = ii;
 		}
 	      retval.maybe_compress ();
-	      // Calculation of 1-norm of inv(*this)
+	      if (calc_cond)
-	      for (octave_idx_type i = 0; i < nm; i++)
+		{
-		work[i] = 0.;
+		  // Calculation of 1-norm of inv(*this)
+		  for (octave_idx_type i = 0; i < nm; i++)
-	      for (octave_idx_type j = 0; j < nr; j++)
+		    work[i] = 0.;
-		{
-		  work[j] = 1.;
+		  for (octave_idx_type j = 0; j < nr; j++)
+		    {
-		  for (octave_idx_type k = j; k < nc; k++)
+		      work[j] = 1.;
-		    {
+		      for (octave_idx_type k = j; k < nc; k++)
-		      if (work[k] != 0.)
 			{
-			  Complex tmp = work[k] / data(cidx(k));
-			  work[k] = tmp;
+			  if (work[k] != 0.)
-			  for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
 			    {
-			      octave_idx_type iidx = ridx(i);
+			      Complex tmp = work[k] / data(cidx(k));
-			      work[iidx] = work[iidx] - tmp * data(i);
+			      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;
-		  double atmp = 0;
+		      for (octave_idx_type i = j; i < nc; i++)
-		  for (octave_idx_type i = j; i < nc; i++)
+			{
-		    {
+			  atmp += std::abs(work[i]);
-		      atmp += std::abs(work[i]);
+			  work[i] = 0.;
-		      work[i] = 0.;
+			}
-		    }
+		      if (atmp > ainvnorm)
-		  if (atmp > ainvnorm)
+			ainvnorm = atmp;
-		    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",
 		   rcond);
 	    }
 	  if (rcond_plus_one == 1.0 || xisnan (rcond))
 	    {
 	      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",
 		   rcond);
 	    }
 return retval;
 }
 ComplexMatrix
-SparseComplexMatrix::trisolve (SparseType &mattype, const Matrix& b, octave_idx_type& err,
+SparseComplexMatrix::trisolve (SparseType &mattype, const Matrix& b,
-			       double& rcond,
+			       octave_idx_type& err, double& rcond,
-			       solve_singularity_handler sing_handler) const
+			       solve_singularity_handler sing_handler,
+			       bool calc_cond) const
 {
 ComplexMatrix retval;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 err = 0;
 if (nr == 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
 volatile int typ = mattype.type ();
 mattype.info ();
 	    {
 	      rcond = 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");
 	    }
 return retval;
 }
 SparseComplexMatrix
 SparseComplexMatrix::trisolve (SparseType &mattype, const SparseMatrix& b,
-		     octave_idx_type& err, double& rcond,
+			       octave_idx_type& err, double& rcond,
-		     solve_singularity_handler sing_handler) const
+			       solve_singularity_handler sing_handler,
+			       bool calc_cond) const
 {
 SparseComplexMatrix retval;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 err = 0;
 if (nr == 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
 int typ = mattype.type ();
 mattype.info ();
 	  if (f77_exception_encountered)
 	    (*current_liboctave_error_handler)
 	      ("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");
 		}
 	      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);
 		  for (volatile octave_idx_type j = 0; j < b_nc; j++)
 		    {
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 err = 0;
 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
 volatile int typ = mattype.type ();
 mattype.info ();
 	    {
 	      rcond = 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");
 	    }
 	}
 return retval;
 }
 SparseComplexMatrix
 SparseComplexMatrix::trisolve (SparseType &mattype,
-		     const SparseComplexMatrix& b, octave_idx_type& err, double& rcond,
+			       const SparseComplexMatrix& b,
-		     solve_singularity_handler sing_handler) const
+			       octave_idx_type& err, double& rcond,
+			       solve_singularity_handler sing_handler,
+			       bool calc_cond) const
 {
 SparseComplexMatrix retval;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 err = 0;
 if (nr == 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
 int typ = mattype.type ();
 mattype.info ();
 	  if (f77_exception_encountered)
 	    (*current_liboctave_error_handler)
 	      ("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");
 		}
 	      else
 		  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 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++)
 return retval;
 }
 ComplexMatrix
-SparseComplexMatrix::bsolve (SparseType &mattype, const Matrix& b, octave_idx_type& err,
+SparseComplexMatrix::bsolve (SparseType &mattype, const Matrix& b,
-			     double& rcond,
+			     octave_idx_type& err, double& rcond,
-			     solve_singularity_handler sing_handler) const
+			     solve_singularity_handler sing_handler,
+			     bool calc_cond) const
 {
 ComplexMatrix retval;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 		if (ri >= j)
 		  m_band(ri - j, j) = data(i);
 	      }
 	  // 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),
 				     nr, n_lower, tmp_data, ldm, err
 				     F77_CHAR_ARG_LEN (1)));
 	  if (f77_exception_encountered)
 	    (*current_liboctave_error_handler)
 	      ("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 ();
 		  typ = SparseType::Banded;
 		  err = 0;
 		}
 	      else
 		{
-		  // Unfortunately, the time to calculate the condition
+		  if (calc_cond)
-		  // number is dominant for narrow banded matrices and
+		    {
-		  // so we rely on the "err" flag from xPBTRF to flag
+		      Array<Complex> z (2 * nr);
-		  // singularity. The commented code below is left here
+		      Complex *pz = z.fortran_vec ();
-		  // for reference
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
-		  //Array<double> z (3 * nr);
-		  //Complex *pz = z.fortran_vec ();
+		      F77_XFCN (zpbcon, ZPBCON,
-		  //Array<octave_idx_type> iz (nr);
+		      	(F77_CONST_CHAR_ARG2 (&job, 1),
-		  //octave_idx_type *piz = iz.fortran_vec ();
+		      	 nr, n_lower, tmp_data, ldm,
-		  //
+		      	 anorm, rcond, pz, piz, err
-		  //F77_XFCN (zpbcon, ZGBCON,
+		      	 F77_CHAR_ARG_LEN (1)));
-		  //	(F77_CONST_CHAR_ARG2 (&job, 1),
-		  //	 nr, n_lower, tmp_data, ldm,
+		      if (f77_exception_encountered)
-		  //	 anorm, rcond, pz, piz, err
+		      	(*current_liboctave_error_handler)
-		  //	 F77_CHAR_ARG_LEN (1)));
+		      	  ("unrecoverable error in zpbcon");
-		  //
-		  //
+		      if (err != 0)
-		  //if (f77_exception_encountered)
+		      	err = -2;
-		  //	(*current_liboctave_error_handler)
-		  //	  ("unrecoverable error in zpbcon");
+		      volatile double rcond_plus_one = rcond + 1.0;
-		  //
-		  //if (err != 0)
+		      if (rcond_plus_one == 1.0 || xisnan (rcond))
-		  //	err = -2;
+		        {
-		  //
+		          err = -2;
-		  //volatile double rcond_plus_one = rcond + 1.0;
-		  //
+		          if (sing_handler)
-		  //if (rcond_plus_one == 1.0 || xisnan (rcond))
+			    {
-		  //  {
+			      sing_handler (rcond);
-		  //    err = -2;
+			      mattype.mark_as_rectangular ();
-		  //
+			    }
-		  //    if (sing_handler)
+		          else
-		  //      sing_handler (rcond);
+		            (*current_liboctave_error_handler)
-		  //    else
+		              ("matrix singular to machine precision, rcond = %g",
-		  //      (*current_liboctave_error_handler)
+		               rcond);
-		  //        ("matrix singular to machine precision, rcond = %g",
+		        }
-		  //         rcond);
+		    }
-		  //  }
+		  else
-		  //else
+		    rcond = 1.0;
-		  //    REST OF CODE, EXCEPT rcond=1
+		  if (err == 0)
-		  rcond = 1.;
+		    {
-		  retval = ComplexMatrix (b);
+		      retval = ComplexMatrix (b);
-		  Complex *result = retval.fortran_vec ();
+		      Complex *result = retval.fortran_vec ();
-		  octave_idx_type b_nc = b.cols ();
+		      octave_idx_type b_nc = b.cols ();
-		  F77_XFCN (zpbtrs, ZPBTRS,
+		      F77_XFCN (zpbtrs, ZPBTRS,
-			    (F77_CONST_CHAR_ARG2 (&job, 1),
+				(F77_CONST_CHAR_ARG2 (&job, 1),
-			     nr, n_lower, b_nc, tmp_data,
+				 nr, n_lower, b_nc, tmp_data,
-			     ldm, result, b.rows(), err
+				 ldm, result, b.rows(), err
-			     F77_CHAR_ARG_LEN (1)));
+				 F77_CHAR_ARG_LEN (1)));
-		  if (f77_exception_encountered)
+		      if (f77_exception_encountered)
-		    (*current_liboctave_error_handler)
+			(*current_liboctave_error_handler)
-		      ("unrecoverable error in zpbtrs");
+			  ("unrecoverable error in zpbtrs");
-		  if (err != 0)
+		      if (err != 0)
-		    {
+			{
-		      (*current_liboctave_error_handler)
+			  (*current_liboctave_error_handler)
-			("SparseMatrix::solve solve failed");
+			    ("SparseMatrix::solve solve failed");
-		      err = -1;
+			  err = -1;
+			}
 		    }
 		}
 	    }
 	}
 	    }
 	  for (octave_idx_type j = 0; j < nc; j++)
 	    for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
 	      m_band(ridx(i) - j + n_lower + n_upper, j) = data(i);
+	  // Calculate the norm of the matrix, for later use.
+	  double anorm;
+	  if (calc_cond)
+	    {
+	      for (octave_idx_type j = 0; j < nr; j++)
+		{
+		  double atmp = 0.;
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		    atmp += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
+	    }
 	  Array<octave_idx_type> ipvt (nr);
 	  octave_idx_type *pipvt = ipvt.fortran_vec ();
 	  F77_XFCN (zgbtrf, ZGBTRF, (nr, nc, n_lower, n_upper, tmp_data,
 	      ("unrecoverable error in zgbtrf");
 	  else
 	    {
 	      // 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';
+		  if (calc_cond)
+		    {
-		  // Unfortunately, the time to calculate the condition
+		      char job = '1';
-		  // number is dominant for narrow banded matrices and
+		      Array<Complex> z (2 * nr);
-		  // so we rely on the "err" flag from xPBTRF to flag
+		      Complex *pz = z.fortran_vec ();
-		  // singularity. The commented code below is left here
+		      Array<double> iz (nr);
-		  // for reference
+		      double *piz = iz.fortran_vec ();
-		  //F77_XFCN (zgbcon, ZGBCON,
+		      F77_XFCN (zgbcon, ZGBCON,
-		  //	(F77_CONST_CHAR_ARG2 (&job, 1),
+		      	(F77_CONST_CHAR_ARG2 (&job, 1),
-		  //	 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
+		      	 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
-		  //	 anorm, rcond, pz, piz, err
+		      	 anorm, rcond, pz, piz, err
-		  //	 F77_CHAR_ARG_LEN (1)));
+		      	 F77_CHAR_ARG_LEN (1)));
-		  //
-		  //if (f77_exception_encountered)
+		      if (f77_exception_encountered)
-		  //  (*current_liboctave_error_handler)
+		        (*current_liboctave_error_handler)
-		  //    ("unrecoverable error in zgbcon");
+		          ("unrecoverable error in zgbcon");
-		  //
-		  // if (err != 0)
+		       if (err != 0)
-		  //  err = -2;
+		        err = -2;
-		  //
-		  //volatile double rcond_plus_one = rcond + 1.0;
+		      volatile double rcond_plus_one = rcond + 1.0;
-		  //
-		  //if (rcond_plus_one == 1.0 || xisnan (rcond))
+		      if (rcond_plus_one == 1.0 || xisnan (rcond))
-		  //  {
+		        {
-		  //    err = -2;
+		          err = -2;
-		  //
-		  //    if (sing_handler)
+		          if (sing_handler)
-		  //      sing_handler (rcond);
+			    {
-		  //    else
+			      sing_handler (rcond);
-		  //      (*current_liboctave_error_handler)
+			      mattype.mark_as_rectangular ();
-		  //        ("matrix singular to machine precision, rcond = %g",
+			    }
-		  //         rcond);
+		          else
-		  //  }
+		            (*current_liboctave_error_handler)
-		  //else
+		              ("matrix singular to machine precision, rcond = %g",
-		  //  REST OF CODE, EXCEPT rcond=1
+		               rcond);
+		        }
-		  rcond = 1.;
+		    }
-		  retval = ComplexMatrix (b);
+		  else
-		  Complex *result = retval.fortran_vec ();
+		    rcond = 1.;
-		  octave_idx_type b_nc = b.cols ();
+		  if (err == 0)
+		    {
-		  job = 'N';
+		      retval = ComplexMatrix (b);
-		  F77_XFCN (zgbtrs, ZGBTRS,
+		      Complex *result = retval.fortran_vec ();
-			    (F77_CONST_CHAR_ARG2 (&job, 1),
-			     nr, n_lower, n_upper, b_nc, tmp_data,
+		      octave_idx_type b_nc = b.cols ();
-			     ldm, pipvt, result, b.rows(), err
-			     F77_CHAR_ARG_LEN (1)));
+		      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)
+		      if (f77_exception_encountered)
-		    (*current_liboctave_error_handler)
+			(*current_liboctave_error_handler)
-		      ("unrecoverable error in zgbtrs");
+			  ("unrecoverable error in zgbtrs");
+		    }
 		}
 	    }
 	}
 else if (typ != SparseType::Banded_Hermitian)
 	(*current_liboctave_error_handler) ("incorrect matrix type");
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	      {
 		octave_idx_type ri = ridx (i);
 		if (ri >= j)
 		  m_band(ri - j, j) = data(i);
 	      }
+	  // Calculate the norm of the matrix, for later use.
+	  double anorm;
+	  if (calc_cond)
+	    anorm = m_band.abs().sum().row(0).max();
 	  char job = 'L';
 	  F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
 				     nr, n_lower, tmp_data, ldm, err
 				     F77_CHAR_ARG_LEN (1)));
 	  if (f77_exception_encountered)
 	    (*current_liboctave_error_handler)
 	      ("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.;
+		  if (calc_cond)
-		  octave_idx_type b_nr = b.rows ();
+		    {
-		  octave_idx_type b_nc = b.cols ();
+		      Array<Complex> z (2 * nr);
-		  OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
+		      Complex *pz = z.fortran_vec ();
+		      Array<double> iz (nr);
-		  // Take a first guess that the number of non-zero terms
+		      double *piz = iz.fortran_vec ();
-		  // will be as many as in b
-		  volatile octave_idx_type x_nz = b.nzmax ();
+		      F77_XFCN (zpbcon, ZPBCON,
-		  volatile octave_idx_type ii = 0;
+		      	(F77_CONST_CHAR_ARG2 (&job, 1),
-		  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
+		      	 nr, n_lower, tmp_data, ldm,
+		      	 anorm, rcond, pz, piz, err
-		  retval.xcidx(0) = 0;
+		      	 F77_CHAR_ARG_LEN (1)));
-		  for (volatile octave_idx_type j = 0; j < b_nc; j++)
-		    {
+		      if (f77_exception_encountered)
-		      for (octave_idx_type i = 0; i < b_nr; i++)
+		      	(*current_liboctave_error_handler)
-			Bx[i] = b.elem (i, j);
+		      	  ("unrecoverable error in zpbcon");
-		      F77_XFCN (zpbtrs, ZPBTRS,
+		      if (err != 0)
-				(F77_CONST_CHAR_ARG2 (&job, 1),
+		      	err = -2;
-				 nr, n_lower, 1, tmp_data,
-				 ldm, Bx, b_nr, err
+		      volatile double rcond_plus_one = rcond + 1.0;
-				 F77_CHAR_ARG_LEN (1)));
+		      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++)
+			{
+			  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)
+			  if (f77_exception_encountered)
-			{
+			    {
-			  (*current_liboctave_error_handler)
+			      (*current_liboctave_error_handler)
-			    ("unrecoverable error in dpbtrs");
+				("unrecoverable error in dpbtrs");
-			  err = -1;
+			      err = -1;
-			  break;
+			      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)
+		      retval.maybe_compress ();
-			{
+		    }
-			  (*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 ();
 		}
 	    }
 	}
 if (typ == SparseType::Banded)
 	    }
 	  for (octave_idx_type j = 0; j < nc; j++)
 	    for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
 	      m_band(ridx(i) - j + n_lower + n_upper, j) = data(i);
+	  // Calculate the norm of the matrix, for later use.
+	  double anorm;
+	  if (calc_cond)
+	    {
+	      for (octave_idx_type j = 0; j < nr; j++)
+		{
+		  double atmp = 0.;
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		    atmp += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
+	    }
 	  Array<octave_idx_type> ipvt (nr);
 	  octave_idx_type *pipvt = ipvt.fortran_vec ();
 	  F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
 	  if (f77_exception_encountered)
 	    (*current_liboctave_error_handler)
 	      ("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");
 		}
 	      else
 		{
-		  char job = 'N';
+		  if (calc_cond)
-		  volatile octave_idx_type x_nz = b.nzmax ();
+		    {
-		  octave_idx_type b_nc = b.cols ();
+		      char job = '1';
-		  retval = SparseComplexMatrix (nr, b_nc, x_nz);
+		      Array<Complex> z (2 * nr);
-		  retval.xcidx(0) = 0;
+		      Complex *pz = z.fortran_vec ();
-		  volatile octave_idx_type ii = 0;
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
-		  OCTAVE_LOCAL_BUFFER (Complex, work, nr);
+		      F77_XFCN (zgbcon, ZGBCON,
-		  for (volatile octave_idx_type j = 0; j < b_nc; j++)
+		      	(F77_CONST_CHAR_ARG2 (&job, 1),
-		    {
+		      	 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
-		      for (octave_idx_type i = 0; i < nr; i++)
+		      	 anorm, rcond, pz, piz, err
-			work[i] = 0.;
+		      	 F77_CHAR_ARG_LEN (1)));
-		      for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
-			work[b.ridx(i)] = b.data(i);
+		      if (f77_exception_encountered)
+		        (*current_liboctave_error_handler)
-		      F77_XFCN (zgbtrs, ZGBTRS,
+		          ("unrecoverable error in zgbcon");
-				(F77_CONST_CHAR_ARG2 (&job, 1),
-				 nr, n_lower, n_upper, 1, tmp_data,
+		       if (err != 0)
-				 ldm, pipvt, work, b.rows (), err
+		        err = -2;
-				 F77_CHAR_ARG_LEN (1)));
+		      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++)
+			{
+			  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)
+			  if (f77_exception_encountered)
-			{
+			    {
-			  (*current_liboctave_error_handler)
+			      (*current_liboctave_error_handler)
-			    ("unrecoverable error in zgbtrs");
+				("unrecoverable error in zgbtrs");
-			  break;
+			      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
+		      retval.maybe_compress ();
-		      // space in retval if needed
+		    }
-		      octave_idx_type new_nnz = 0;
-		      for (octave_idx_type i = 0; i < nr; i++)
-			if (work[i] != 0.)
-			  new_nnz++;
-		      if (ii + new_nnz > x_nz)
-			{
-			  // Resize the sparse matrix
-			  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
-			  retval.change_capacity (sz);
-			  x_nz = sz;
-			}
-		      for (octave_idx_type i = 0; i < nr; i++)
-			if (work[i] != 0.)
-			  {
-			    retval.xridx(ii) = i;
-			    retval.xdata(ii++) = work[i];
-			  }
-		      retval.xcidx(j+1) = ii;
-		    }
-		  retval.maybe_compress ();
 		}
 	    }
 	}
 else if (typ != SparseType::Banded_Hermitian)
 	(*current_liboctave_error_handler) ("incorrect matrix type");
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	      {
 		octave_idx_type ri = ridx (i);
 		if (ri >= j)
 		  m_band(ri - j, j) = data(i);
 	      }
+	  // Calculate the norm of the matrix, for later use.
+	  double anorm;
+	  if (calc_cond)
+	    anorm = m_band.abs().sum().row(0).max();
 	  char job = 'L';
 	  F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
 				     nr, n_lower, tmp_data, ldm, err
 				     F77_CHAR_ARG_LEN (1)));
 	  if (f77_exception_encountered)
 	    (*current_liboctave_error_handler)
 	      ("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.;
+		  if (calc_cond)
-		  octave_idx_type b_nr = b.rows ();
+		    {
-		  octave_idx_type b_nc = b.cols ();
+		      Array<Complex> z (2 * nr);
-		  retval = ComplexMatrix (b);
+		      Complex *pz = z.fortran_vec ();
-		  Complex *result = retval.fortran_vec ();
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
-		  F77_XFCN (zpbtrs, ZPBTRS,
-			    (F77_CONST_CHAR_ARG2 (&job, 1),
+		      F77_XFCN (zpbcon, ZPBCON,
-			     nr, n_lower, b_nc, tmp_data,
+		      	(F77_CONST_CHAR_ARG2 (&job, 1),
-			     ldm, result, b_nr, err
+		      	 nr, n_lower, tmp_data, ldm,
-			     F77_CHAR_ARG_LEN (1)));
+		      	 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)
+		      if (f77_exception_encountered)
-		    {
+			{
-		      (*current_liboctave_error_handler)
+			  (*current_liboctave_error_handler)
-			("unrecoverable error in zpbtrs");
+			    ("unrecoverable error in zpbtrs");
-		      err = -1;
+			  err = -1;
-		    }
+			}
-		  if (err != 0)
+		      if (err != 0)
-		    {
+			{
-		      (*current_liboctave_error_handler)
+			  (*current_liboctave_error_handler)
-			("SparseComplexMatrix::solve solve failed");
+			    ("SparseComplexMatrix::solve solve failed");
-		      err = -1;
+			  err = -1;
+			}
 		    }
 		}
 	    }
 	}
 	    }
 	  for (octave_idx_type j = 0; j < nc; j++)
 	    for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
 	      m_band(ridx(i) - j + n_lower + n_upper, j) = data(i);
+	  // Calculate the norm of the matrix, for later use.
+	  double anorm;
+	  if (calc_cond)
+	    {
+	      for (octave_idx_type j = 0; j < nr; j++)
+		{
+		  double atmp = 0.;
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		    atmp += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
+	    }
 	  Array<octave_idx_type> ipvt (nr);
 	  octave_idx_type *pipvt = ipvt.fortran_vec ();
 	  F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
 	  if (f77_exception_encountered)
 	    (*current_liboctave_error_handler)
 	      ("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';
+		  if (calc_cond)
-		  octave_idx_type b_nc = b.cols ();
+		    {
-		  retval = ComplexMatrix (b);
+		      char job = '1';
-		  Complex *result = retval.fortran_vec ();
+		      Array<Complex> z (2 * nr);
+		      Complex *pz = z.fortran_vec ();
-		  F77_XFCN (zgbtrs, ZGBTRS,
+		      Array<double> iz (nr);
-			    (F77_CONST_CHAR_ARG2 (&job, 1),
+		      double *piz = iz.fortran_vec ();
-			     nr, n_lower, n_upper, b_nc, tmp_data,
-			     ldm, pipvt, result, b.rows (), err
+		      F77_XFCN (zgbcon, ZGBCON,
-			     F77_CHAR_ARG_LEN (1)));
+		      	(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)
+		      if (f77_exception_encountered)
-		    {
+			{
-		      (*current_liboctave_error_handler)
+			  (*current_liboctave_error_handler)
-			("unrecoverable error in dgbtrs");
+			    ("unrecoverable error in dgbtrs");
+			}
 		    }
 		}
 	    }
 	}
 else if (typ != SparseType::Banded_Hermitian)
 return retval;
 }
 SparseComplexMatrix
 SparseComplexMatrix::bsolve (SparseType &mattype, const SparseComplexMatrix& b,
-		     octave_idx_type& err, double& rcond,
+			     octave_idx_type& err, double& rcond,
-		     solve_singularity_handler sing_handler) const
+			     solve_singularity_handler sing_handler,
+			     bool calc_cond) const
 {
 SparseComplexMatrix retval;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	      {
 		octave_idx_type ri = ridx (i);
 		if (ri >= j)
 		  m_band(ri - j, j) = data(i);
 	      }
+	  // Calculate the norm of the matrix, for later use.
+	  double anorm;
+	  if (calc_cond)
+	    anorm = m_band.abs().sum().row(0).max();
 	  char job = 'L';
 	  F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
 				     nr, n_lower, tmp_data, ldm, err
 				     F77_CHAR_ARG_LEN (1)));
 	  if (f77_exception_encountered)
 	    (*current_liboctave_error_handler)
 	      ("unrecoverable error in zpbtrf");
 	  else
 	    {
-	      rcond = 0.0;
 	      if (err != 0)
 		{
 		  // Matrix is not positive definite!! Fall through to
 		  // unsymmetric banded solver.
 		  mattype.mark_as_unsymmetric ();
 		  typ = SparseType::Banded;
+		  rcond = 0.0;
 		  err = 0;
 		}
 	      else
 		{
-		  rcond = 1.;
+		  if (calc_cond)
-		  octave_idx_type b_nr = b.rows ();
+		    {
-		  octave_idx_type b_nc = b.cols ();
+		      Array<Complex> z (2 * nr);
-		  OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
+		      Complex *pz = z.fortran_vec ();
+		      Array<double> iz (nr);
-		  // Take a first guess that the number of non-zero terms
+		      double *piz = iz.fortran_vec ();
-		  // will be as many as in b
-		  volatile octave_idx_type x_nz = b.nzmax ();
+		      F77_XFCN (zpbcon, ZPBCON,
-		  volatile octave_idx_type ii = 0;
+		      	(F77_CONST_CHAR_ARG2 (&job, 1),
-		  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
+		      	 nr, n_lower, tmp_data, ldm,
+		      	 anorm, rcond, pz, piz, err
-		  retval.xcidx(0) = 0;
+		      	 F77_CHAR_ARG_LEN (1)));
-		  for (volatile octave_idx_type j = 0; j < b_nc; j++)
-		    {
+		      if (f77_exception_encountered)
+		      	(*current_liboctave_error_handler)
-		      for (octave_idx_type i = 0; i < b_nr; i++)
+		      	  ("unrecoverable error in zpbcon");
-			Bx[i] = b (i,j);
+		      if (err != 0)
-		      F77_XFCN (zpbtrs, ZPBTRS,
+		      	err = -2;
-				(F77_CONST_CHAR_ARG2 (&job, 1),
-				 nr, n_lower, 1, tmp_data,
+		      volatile double rcond_plus_one = rcond + 1.0;
-				 ldm, Bx, b_nr, err
-				 F77_CHAR_ARG_LEN (1)));
+		      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++)
+			{
+			  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)
+			  if (f77_exception_encountered)
-			{
+			    {
-			  (*current_liboctave_error_handler)
+			      (*current_liboctave_error_handler)
-			    ("unrecoverable error in zpbtrs");
+				("unrecoverable error in zpbtrs");
-			  err = -1;
+			      err = -1;
-			  break;
+			      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;
+			    }
+			  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;
 			}
-		      if (err != 0)
+		      retval.maybe_compress ();
-			{
+		    }
-			  (*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;
-			}
-		      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 ();
 		}
 	    }
 	}
 if (typ == SparseType::Banded)
 	    }
 	  for (octave_idx_type j = 0; j < nc; j++)
 	    for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
 	      m_band(ridx(i) - j + n_lower + n_upper, j) = data(i);
+	  // Calculate the norm of the matrix, for later use.
+	  double anorm;
+	  if (calc_cond)
+	    {
+	      for (octave_idx_type j = 0; j < nr; j++)
+		{
+		  double atmp = 0.;
+		  for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
+		    atmp += std::abs(data(i));
+		  if (atmp > anorm)
+		    anorm = atmp;
+		}
+	    }
 	  Array<octave_idx_type> ipvt (nr);
 	  octave_idx_type *pipvt = ipvt.fortran_vec ();
 	  F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
 	  if (f77_exception_encountered)
 	    (*current_liboctave_error_handler)
 	      ("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");
 		}
 	      else
 		{
-		  char job = 'N';
+		  if (calc_cond)
-		  volatile octave_idx_type x_nz = b.nzmax ();
+		    {
-		  octave_idx_type b_nc = b.cols ();
+		      char job = '1';
-		  retval = SparseComplexMatrix (nr, b_nc, x_nz);
+		      Array<Complex> z (2 * nr);
-		  retval.xcidx(0) = 0;
+		      Complex *pz = z.fortran_vec ();
-		  volatile octave_idx_type ii = 0;
+		      Array<double> iz (nr);
+		      double *piz = iz.fortran_vec ();
-		  OCTAVE_LOCAL_BUFFER (Complex, Bx, nr);
+		      F77_XFCN (zgbcon, ZGBCON,
-		  for (volatile octave_idx_type j = 0; j < b_nc; j++)
+		      	(F77_CONST_CHAR_ARG2 (&job, 1),
-		    {
+		      	 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
-		      for (octave_idx_type i = 0; i < nr; i++)
+		      	 anorm, rcond, pz, piz, err
-			Bx[i] = 0.;
+		      	 F77_CHAR_ARG_LEN (1)));
-		      for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
+		      if (f77_exception_encountered)
-			Bx[b.ridx(i)] = b.data(i);
+		        (*current_liboctave_error_handler)
+		          ("unrecoverable error in zgbcon");
-		      F77_XFCN (zgbtrs, ZGBTRS,
-				(F77_CONST_CHAR_ARG2 (&job, 1),
+		       if (err != 0)
-				 nr, n_lower, n_upper, 1, tmp_data,
+		        err = -2;
-				 ldm, pipvt, Bx, b.rows (), err
-				 F77_CHAR_ARG_LEN (1)));
+		      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++)
+			{
+			  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)
+			  if (f77_exception_encountered)
-			{
+			    {
-			  (*current_liboctave_error_handler)
+			      (*current_liboctave_error_handler)
-			    ("unrecoverable error in dgbtrs");
+				("unrecoverable error in dgbtrs");
-			  break;
+			      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
+		      retval.maybe_compress ();
-		      // 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 ();
 		}
 	    }
 	}
 else if (typ != SparseType::Banded_Hermitian)
 	(*current_liboctave_error_handler) ("incorrect matrix type");
 return retval;
 }
 void *
-SparseComplexMatrix::factorize (octave_idx_type& err, double &rcond, Matrix &Control,
+SparseComplexMatrix::factorize (octave_idx_type& err, double &rcond,
-				Matrix &Info,
+				Matrix &Control, Matrix &Info,
-				solve_singularity_handler sing_handler) const
+				solve_singularity_handler sing_handler,
+				bool calc_cond) const
 {
 // The return values
 void *Numeric = 0;
 err = 0;
 status = UMFPACK_ZNAME (numeric) (Ap, Ai,
 				   X_CAST (const double *, Ax), NULL,
 				   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 ||
 	  rcond_plus_one == 1.0 || xisnan (rcond))
 	{
 return Numeric;
 }
 ComplexMatrix
-SparseComplexMatrix::fsolve (SparseType &mattype, const Matrix& b, octave_idx_type& err,
+SparseComplexMatrix::fsolve (SparseType &mattype, const Matrix& b,
-			     double& rcond,
+			     octave_idx_type& err, double& rcond,
-			     solve_singularity_handler sing_handler) const
+			     solve_singularity_handler sing_handler,
+			     bool calc_cond) const
 {
 ComplexMatrix retval;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	  cholmod_factor *L;
 	  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)
 	    {
 	      // Either its indefinite or singular. Try UMFPACK
 	      if (rcond_plus_one == 1.0 || xisnan (rcond))
 		{
 		  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",
 		       rcond);
 if (typ == SparseType::Full)
 	{
 #ifdef HAVE_UMFPACK
 	  Matrix Control, Info;
 	  void *Numeric = factorize (err, rcond, Control, Info,
-				     sing_handler);
+				     sing_handler, calc_cond);
 	  if (err == 0)
 	    {
 	      octave_idx_type b_nr = b.rows ();
 	      octave_idx_type b_nc = b.cols ();
 	      UMFPACK_ZNAME (report_info) (control, info);
 	      UMFPACK_ZNAME (free_numeric) (&Numeric);
 	    }
+	  else
+	    mattype.mark_as_rectangular ();
 #else
 	  (*current_liboctave_error_handler) ("UMFPACK not installed");
 #endif
 	}
 else if (typ != SparseType::Hermitian)
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	  cholmod_factor *L;
 	  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)
 	    {
 	      // Either its indefinite or singular. Try UMFPACK
 	      if (rcond_plus_one == 1.0 || xisnan (rcond))
 		{
 		  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",
 		       rcond);
 if (typ == SparseType::Full)
 	{
 #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)
 	    {
 	      octave_idx_type b_nr = b.rows ();
 	      octave_idx_type b_nc = b.cols ();
 	      OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr);
 #endif
 	      // 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);
 	      OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
 	      UMFPACK_ZNAME (report_info) (control, info);
 	      UMFPACK_ZNAME (free_numeric) (&Numeric);
 	    }
+	  else
+	    mattype.mark_as_rectangular ();
 #else
 	  (*current_liboctave_error_handler) ("UMFPACK not installed");
 #endif
 	}
 else if (typ != SparseType::Hermitian)
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	  cholmod_factor *L;
 	  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)
 	    {
 	      // Either its indefinite or singular. Try UMFPACK
 	      if (rcond_plus_one == 1.0 || xisnan (rcond))
 		{
 		  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",
 		       rcond);
 if (typ == SparseType::Full)
 	{
 #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)
 	    {
 	      octave_idx_type b_nr = b.rows ();
 	      octave_idx_type b_nc = b.cols ();
 	      UMFPACK_ZNAME (report_info) (control, info);
 	      UMFPACK_ZNAME (free_numeric) (&Numeric);
 	    }
+	  else
+	    mattype.mark_as_rectangular ();
 #else
 	  (*current_liboctave_error_handler) ("UMFPACK not installed");
 #endif
 	}
 else if (typ != SparseType::Hermitian)
 }
 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;
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 	  cholmod_factor *L;
 	  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)
 	    {
 	      // Either its indefinite or singular. Try UMFPACK
 	      if (rcond_plus_one == 1.0 || xisnan (rcond))
 		{
 		  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",
 		       rcond);
 if (typ == SparseType::Full)
 	{
 #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)
 	    {
 	      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
-	      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);
 	      OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
 	      UMFPACK_ZNAME (report_info) (control, info);
 	      UMFPACK_ZNAME (free_numeric) (&Numeric);
 	    }
+	  else
+	    mattype.mark_as_rectangular ();
 #else
 	  (*current_liboctave_error_handler) ("UMFPACK not installed");
 #endif
 	}
 else if (typ != SparseType::Hermitian)
 ComplexMatrix
 SparseComplexMatrix::solve (SparseType &mattype, const Matrix& b, 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);
+retval = fsolve (mattype, b, err, rcond, sing_handler, true);
-else
+else if (typ != SparseType::Rectangular)
 {
-(*current_liboctave_error_handler)
+(*current_liboctave_error_handler) ("unknown matrix type");
-	("matrix dimension mismatch solution of linear equations");
 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
 SparseComplexMatrix::solve (SparseType &mattype, const SparseMatrix& b) const
 {
 SparseComplexMatrix
 SparseComplexMatrix::solve (SparseType &mattype, const SparseMatrix& b,
 			    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);
+retval = fsolve (mattype, b, err, rcond, sing_handler, true);
-else
+else if (typ != SparseType::Rectangular)
 {
-(*current_liboctave_error_handler)
+(*current_liboctave_error_handler) ("unknown matrix type");
-	("matrix dimension mismatch solution of linear equations");
 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
 SparseComplexMatrix::solve (SparseType &mattype, const ComplexMatrix& b) const
 {
 ComplexMatrix
 SparseComplexMatrix::solve (SparseType &mattype, const ComplexMatrix& b,
 		     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);
+retval = fsolve (mattype, b, err, rcond, sing_handler, true);
-else
+else if (typ != SparseType::Rectangular)
 {
-(*current_liboctave_error_handler)
+(*current_liboctave_error_handler) ("unknown matrix type");
-	("matrix dimension mismatch solution of linear equations");
 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
 SparseComplexMatrix::solve (SparseType &mattype,
 			    const SparseComplexMatrix& b) const
 SparseComplexMatrix
 SparseComplexMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b,
 			    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);
+retval = fsolve (mattype, b, err, rcond, sing_handler, true);
-else
+else if (typ != SparseType::Rectangular)
 {
-(*current_liboctave_error_handler)
+(*current_liboctave_error_handler) ("unknown matrix type");
-	("matrix dimension mismatch solution of linear equations");
 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
 SparseComplexMatrix::solve (SparseType &mattype, const ColumnVector& b) const
 {
 {
 ComplexMatrix tmp (b);
 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);
 octave_idx_type ii = 0;
 SparseComplexMatrix
 SparseComplexMatrix::map (c_c_Mapper f) const
 {
 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
 octave_idx_type nel = (f_zero ? 0 : nr*nc - nz);
 for (octave_idx_type i = 0; i < nz; i++)
 SparseMatrix
 SparseComplexMatrix::map (d_c_Mapper f) const
 {
 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
 octave_idx_type nel = (f_zero ? 0 : nr*nc - nz);
 for (octave_idx_type i = 0; i < nz; i++)
 SparseBoolMatrix
 SparseComplexMatrix::map (b_c_Mapper f) const
 {
 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
 octave_idx_type nel = (f_zero ? 0 : nr*nc - nz);
 for (octave_idx_type i = 0; i < nz; i++)
 }
 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++)
 {
 Complex val = data (i);
 if (xisinf (val) || xisnan (val))
 // Return true if no elements have imaginary components.
 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++)
 {
 double ip = std::imag (data (i));
 // imaginary) values and return them in MAX_VAL and MIN_VAL.
 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;
 max_val = std::real(data (0));
 }
 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++)
 {
 	Complex val = data (i);
 #undef COL_EXPR
 }
 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);
 for (octave_idx_type i = 0; i < nc + 1; i++)
 }
 SparseComplexMatrix
 operator * (const SparseComplexMatrix& m, const SparseMatrix& a)
 {
-SparseComplexMatrix tmp (a);
+SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, double);
-return m * tmp;
 }
 SparseComplexMatrix
 operator * (const SparseMatrix& m, const SparseComplexMatrix& a)
 {
-SparseComplexMatrix tmp (m);
+SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, Complex);
-return tmp * a;
 }
 SparseComplexMatrix
 operator * (const SparseComplexMatrix& m, const SparseComplexMatrix& a)
 {
-#ifdef HAVE_SPARSE_BLAS
+SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, Complex);
-// XXX FIXME XXX Isn't there a sparse BLAS ??
-#else
-// Use Andy's sparse matrix multiply function
-SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex);
-#endif
 }
 ComplexMatrix
 operator * (const ComplexMatrix& m, const SparseMatrix& a)
 {
-SparseComplexMatrix tmp (a);
+FULL_SPARSE_MUL (ComplexMatrix, double, Complex (0.,0.));
-return m * tmp;
 }
 ComplexMatrix
 operator * (const Matrix& m, const SparseComplexMatrix& a)
 {
-ComplexMatrix tmp (m);
+FULL_SPARSE_MUL (ComplexMatrix, Complex, Complex (0.,0.));
-return tmp * a;
 }
 ComplexMatrix
 operator * (const ComplexMatrix& m, const SparseComplexMatrix& a)
 {
-#ifdef HAVE_SPARSE_BLAS
+FULL_SPARSE_MUL (ComplexMatrix, Complex, Complex (0.,0.));
-// XXX FIXME XXX Isn't there a sparse BLAS ??
-#else
-FULL_SPARSE_MUL (ComplexMatrix, Complex);
-#endif
 }
 ComplexMatrix
 operator * (const SparseComplexMatrix& m, const Matrix& a)
 {
-ComplexMatrix tmp (a);
+SPARSE_FULL_MUL (ComplexMatrix, double, Complex (0.,0.));
-return m * tmp;
 }
 ComplexMatrix
 operator * (const SparseMatrix& m, const ComplexMatrix& a)
 {
-SparseComplexMatrix tmp (m);
+SPARSE_FULL_MUL (ComplexMatrix, Complex, Complex (0.,0.));
-return tmp * a;
 }
 ComplexMatrix
 operator * (const SparseComplexMatrix& m, const ComplexMatrix& a)
 {
-#ifdef HAVE_SPARSE_BLAS
+SPARSE_FULL_MUL (ComplexMatrix, Complex, Complex (0.,0.));
-// XXX FIXME XXX Isn't there a sparse BLAS ??
-#else
-SPARSE_FULL_MUL (ComplexMatrix, Complex);
-#endif
 }
 // XXX FIXME XXX -- it would be nice to share code among the min/max
 // functions below.
 octave_idx_type a_nc = a.cols ();
 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;
 	  for (octave_idx_type i = 0 ; i < a_nc ; i++)
 	    {
 octave_idx_type b_nr = b.rows ();
 octave_idx_type b_nc = b.cols ();
 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;
 	  for (octave_idx_type i = 0 ; i < a_nc ; i++)
 	    {

Mercurial > hg > octave-nkf

comparison liboctave/CSparse.cc @ 5681:233d98d95659