diff liboctave/array/dSparse.cc @ 21136:7cac4e7458f2

maint: clean up code around calls to current_liboctave_error_handler. Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code. * Array-util.cc, Array.cc, CColVector.cc, CDiagMatrix.cc, CMatrix.cc, CNDArray.cc, CRowVector.cc, CSparse.cc, DiagArray2.cc, MArray.cc, PermMatrix.cc, Sparse.cc, Sparse.h, chMatrix.cc, chNDArray.cc, dColVector.cc, dDiagMatrix.cc, dMatrix.cc, dNDArray.cc, dRowVector.cc, dSparse.cc, fCColVector.cc, fCDiagMatrix.cc, fCMatrix.cc, fCNDArray.cc, fCRowVector.cc, fColVector.cc, fDiagMatrix.cc, fMatrix.cc, fNDArray.cc, fRowVector.cc, idx-vector.cc, CmplxAEPBAL.cc, CmplxCHOL.cc, CmplxGEPBAL.cc, CmplxHESS.cc, CmplxLU.cc, CmplxQR.cc, CmplxSCHUR.cc, CmplxSVD.cc, DASPK.cc, EIG.cc, LSODE.cc, Quad.cc, SparseCmplxCHOL.cc, SparseCmplxLU.cc, SparseCmplxQR.cc, SparseQR.cc, SparsedbleCHOL.cc, SparsedbleLU.cc, base-lu.cc, bsxfun-defs.cc, dbleAEPBAL.cc, dbleCHOL.cc, dbleGEPBAL.cc, dbleHESS.cc, dbleLU.cc, dbleQR.cc, dbleSCHUR.cc, dbleSVD.cc, eigs-base.cc, fCmplxAEPBAL.cc, fCmplxCHOL.cc, fCmplxLU.cc, fCmplxQR.cc, fCmplxSCHUR.cc, fEIG.cc, floatAEPBAL.cc, floatCHOL.cc, floatGEPBAL.cc, floatHESS.cc, floatLU.cc, floatQR.cc, floatSCHUR.cc, floatSVD.cc, lo-specfun.cc, oct-fftw.cc, oct-rand.cc, oct-spparms.cc, sparse-base-chol.cc, sparse-dmsolve.cc, file-ops.cc, lo-sysdep.cc, mach-info.cc, oct-env.cc, oct-syscalls.cc, cmd-edit.cc, cmd-hist.cc, data-conv.cc, lo-ieee.cc, lo-regexp.cc, oct-base64.cc, oct-shlib.cc, pathsearch.cc, singleton-cleanup.cc, sparse-util.cc, unwind-prot.cc: Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code.
author Rik <rik@octave.org>
date Sat, 23 Jan 2016 13:52:03 -0800
parents 499b851fbfae
children 623fc7d08cc6
line wrap: on
line diff
--- a/liboctave/array/dSparse.cc	Fri Jan 22 13:45:21 2016 -0500
+++ b/liboctave/array/dSparse.cc	Sat Jan 23 13:52:03 2016 -0800
@@ -751,73 +751,71 @@
 SparseMatrix
 atan2 (const SparseMatrix& x, const SparseMatrix& y)
 {
+  if ((x.rows () != y.rows ()) || (x.cols () != y.cols ()))
+    (*current_liboctave_error_handler) ("matrix size mismatch");
+
+  octave_idx_type x_nr = x.rows ();
+  octave_idx_type x_nc = x.cols ();
+
+  octave_idx_type y_nr = y.rows ();
+  octave_idx_type y_nc = y.cols ();
+
+  if (x_nr != y_nr || x_nc != y_nc)
+    err_nonconformant ("atan2", x_nr, x_nc, y_nr, y_nc);
+
   SparseMatrix r;
 
-  if ((x.rows () == y.rows ()) && (x.cols () == y.cols ()))
+  r = SparseMatrix (x_nr, x_nc, (x.nnz () + y.nnz ()));
+
+  octave_idx_type jx = 0;
+  r.cidx (0) = 0;
+  for (octave_idx_type i = 0 ; i < x_nc ; i++)
     {
-      octave_idx_type x_nr = x.rows ();
-      octave_idx_type x_nc = x.cols ();
-
-      octave_idx_type y_nr = y.rows ();
-      octave_idx_type y_nc = y.cols ();
-
-      if (x_nr != y_nr || x_nc != y_nc)
-        err_nonconformant ("atan2", x_nr, x_nc, y_nr, y_nc);
-
-      r = SparseMatrix (x_nr, x_nc, (x.nnz () + y.nnz ()));
-
-      octave_idx_type jx = 0;
-      r.cidx (0) = 0;
-      for (octave_idx_type i = 0 ; i < x_nc ; i++)
+      octave_idx_type  ja = x.cidx (i);
+      octave_idx_type  ja_max = x.cidx (i+1);
+      bool ja_lt_max= ja < ja_max;
+
+      octave_idx_type  jb = y.cidx (i);
+      octave_idx_type  jb_max = y.cidx (i+1);
+      bool jb_lt_max = jb < jb_max;
+
+      while (ja_lt_max || jb_lt_max)
         {
-          octave_idx_type  ja = x.cidx (i);
-          octave_idx_type  ja_max = x.cidx (i+1);
-          bool ja_lt_max= ja < ja_max;
-
-          octave_idx_type  jb = y.cidx (i);
-          octave_idx_type  jb_max = y.cidx (i+1);
-          bool jb_lt_max = jb < jb_max;
-
-          while (ja_lt_max || jb_lt_max)
-            {
-              octave_quit ();
-              if ((! jb_lt_max)
-                  || (ja_lt_max && (x.ridx (ja) < y.ridx (jb))))
+          octave_quit ();
+          if ((! jb_lt_max)
+              || (ja_lt_max && (x.ridx (ja) < y.ridx (jb))))
+            {
+              r.ridx (jx) = x.ridx (ja);
+              r.data (jx) = atan2 (x.data (ja), 0.);
+              jx++;
+              ja++;
+              ja_lt_max= ja < ja_max;
+            }
+          else if ((! ja_lt_max)
+                   || (jb_lt_max && (y.ridx (jb) < x.ridx (ja))))
+            {
+              jb++;
+              jb_lt_max= jb < jb_max;
+            }
+          else
+            {
+              double tmp = atan2 (x.data (ja), y.data (jb));
+              if (tmp != 0.)
                 {
+                  r.data (jx) = tmp;
                   r.ridx (jx) = x.ridx (ja);
-                  r.data (jx) = atan2 (x.data (ja), 0.);
                   jx++;
-                  ja++;
-                  ja_lt_max= ja < ja_max;
-                }
-              else if ((! ja_lt_max)
-                       || (jb_lt_max && (y.ridx (jb) < x.ridx (ja))))
-                {
-                  jb++;
-                  jb_lt_max= jb < jb_max;
                 }
-              else
-                {
-                  double tmp = atan2 (x.data (ja), y.data (jb));
-                  if (tmp != 0.)
-                    {
-                      r.data (jx) = tmp;
-                      r.ridx (jx) = x.ridx (ja);
-                      jx++;
-                    }
-                  ja++;
-                  ja_lt_max= ja < ja_max;
-                  jb++;
-                  jb_lt_max= jb < jb_max;
-                }
-            }
-          r.cidx (i+1) = jx;
+              ja++;
+              ja_lt_max= ja < ja_max;
+              jb++;
+              jb_lt_max= jb < jb_max;
+            }
         }
-
-      r.maybe_compress ();
+      r.cidx (i+1) = jx;
     }
-  else
-    (*current_liboctave_error_handler) ("matrix size mismatch");
+
+  r.maybe_compress ();
 
   return r;
 }
@@ -859,44 +857,40 @@
 
   if (nr == 0 || nc == 0 || nr != nc)
     (*current_liboctave_error_handler) ("inverse requires square matrix");
+
+  // Print spparms("spumoni") info if requested
+  int typ = mattyp.type ();
+  mattyp.info ();
+
+  if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal)
+    (*current_liboctave_error_handler) ("incorrect matrix type");
+
+  if (typ == MatrixType::Permuted_Diagonal)
+    retval = transpose ();
   else
+    retval = *this;
+
+  // Force make_unique to be called
+  double *v = retval.data ();
+
+  if (calccond)
     {
-      // Print spparms("spumoni") info if requested
-      int typ = mattyp.type ();
-      mattyp.info ();
-
-      if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
+      double dmax = 0.;
+      double dmin = octave_Inf;
+      for (octave_idx_type i = 0; i < nr; i++)
         {
-          if (typ == MatrixType::Permuted_Diagonal)
-            retval = transpose ();
-          else
-            retval = *this;
-
-          // Force make_unique to be called
-          double *v = retval.data ();
-
-          if (calccond)
-            {
-              double dmax = 0.;
-              double dmin = octave_Inf;
-              for (octave_idx_type i = 0; i < nr; i++)
-                {
-                  double tmp = fabs (v[i]);
-                  if (tmp > dmax)
-                    dmax = tmp;
-                  if (tmp < dmin)
-                    dmin = tmp;
-                }
-              rcond = dmin / dmax;
-            }
-
-          for (octave_idx_type i = 0; i < nr; i++)
-            v[i] = 1.0 / v[i];
+          double tmp = fabs (v[i]);
+          if (tmp > dmax)
+            dmax = tmp;
+          if (tmp < dmin)
+            dmin = tmp;
         }
-      else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+      rcond = dmin / dmax;
     }
 
+  for (octave_idx_type i = 0; i < nr; i++)
+    v[i] = 1.0 / v[i];
+
   return retval;
 }
 
@@ -913,262 +907,239 @@
 
   if (nr == 0 || nc == 0 || nr != nc)
     (*current_liboctave_error_handler) ("inverse requires square matrix");
-  else
+
+  // Print spparms("spumoni") info if requested
+  int typ = mattyp.type ();
+  mattyp.info ();
+
+  if (typ != MatrixType::Upper && typ != MatrixType::Permuted_Upper
+      && typ != MatrixType::Lower && typ != MatrixType::Permuted_Lower)
+    (*current_liboctave_error_handler) ("incorrect matrix type");
+
+  double anorm = 0.;
+  double ainvnorm = 0.;
+
+  if (calccond)
     {
-      // Print spparms("spumoni") info if requested
-      int typ = mattyp.type ();
-      mattyp.info ();
-
-      if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper
-          || typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
+      // Calculate the 1-norm of matrix for rcond calculation
+      for (octave_idx_type j = 0; j < nr; j++)
         {
-          double anorm = 0.;
-          double ainvnorm = 0.;
-
-          if (calccond)
-            {
-              // Calculate the 1-norm of matrix for rcond calculation
-              for (octave_idx_type j = 0; j < nr; j++)
-                {
-                  double atmp = 0.;
-                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
-                    atmp += fabs (data (i));
-                  if (atmp > anorm)
-                    anorm = atmp;
-                }
-            }
-
-          if (typ == MatrixType::Upper || typ == MatrixType::Lower)
-            {
-              octave_idx_type nz = nnz ();
-              octave_idx_type cx = 0;
-              octave_idx_type nz2 = nz;
-              retval = SparseMatrix (nr, nc, nz2);
-
-              for (octave_idx_type i = 0; i < nr; i++)
+          double atmp = 0.;
+          for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+            atmp += fabs (data (i));
+          if (atmp > anorm)
+            anorm = atmp;
+        }
+    }
+
+  if (typ == MatrixType::Upper || typ == MatrixType::Lower)
+    {
+      octave_idx_type nz = nnz ();
+      octave_idx_type cx = 0;
+      octave_idx_type nz2 = nz;
+      retval = SparseMatrix (nr, nc, nz2);
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        {
+          octave_quit ();
+          // place the 1 in the identity position
+          octave_idx_type cx_colstart = cx;
+
+          if (cx == nz2)
+            {
+              nz2 *= 2;
+              retval.change_capacity (nz2);
+            }
+
+          retval.xcidx (i) = cx;
+          retval.xridx (cx) = i;
+          retval.xdata (cx) = 1.0;
+          cx++;
+
+          // iterate accross columns of input matrix
+          for (octave_idx_type j = i+1; j < nr; j++)
+            {
+              double v = 0.;
+              // iterate to calculate sum
+              octave_idx_type colXp = retval.xcidx (i);
+              octave_idx_type colUp = cidx (j);
+              octave_idx_type rpX, rpU;
+
+              if (cidx (j) == cidx (j+1))
+                (*current_liboctave_error_handler) ("division by zero");
+
+              do
                 {
                   octave_quit ();
-                  // place the 1 in the identity position
-                  octave_idx_type cx_colstart = cx;
-
+                  rpX = retval.xridx (colXp);
+                  rpU = ridx (colUp);
+
+                  if (rpX < rpU)
+                    colXp++;
+                  else if (rpX > rpU)
+                    colUp++;
+                  else
+                    {
+                      v -= retval.xdata (colXp) * data (colUp);
+                      colXp++;
+                      colUp++;
+                    }
+                }
+              while (rpX < j && rpU < j && colXp < cx && colUp < nz);
+
+              // get A(m,m)
+              if (typ == MatrixType::Upper)
+                colUp = cidx (j+1) - 1;
+              else
+                colUp = cidx (j);
+              double pivot = data (colUp);
+              if (pivot == 0. || ridx (colUp) != j)
+                (*current_liboctave_error_handler) ("division by zero");
+
+              if (v != 0.)
+                {
                   if (cx == nz2)
                     {
                       nz2 *= 2;
                       retval.change_capacity (nz2);
                     }
 
-                  retval.xcidx (i) = cx;
-                  retval.xridx (cx) = i;
-                  retval.xdata (cx) = 1.0;
+                  retval.xridx (cx) = j;
+                  retval.xdata (cx) = v / pivot;
                   cx++;
-
-                  // iterate accross columns of input matrix
-                  for (octave_idx_type j = i+1; j < nr; j++)
-                    {
-                      double v = 0.;
-                      // iterate to calculate sum
-                      octave_idx_type colXp = retval.xcidx (i);
-                      octave_idx_type colUp = cidx (j);
-                      octave_idx_type rpX, rpU;
-
-                      if (cidx (j) == cidx (j+1))
-                        {
-                          (*current_liboctave_error_handler)
-                            ("division by zero");
-                          goto inverse_singular;
-                        }
-
-                      do
-                        {
-                          octave_quit ();
-                          rpX = retval.xridx (colXp);
-                          rpU = ridx (colUp);
-
-                          if (rpX < rpU)
-                            colXp++;
-                          else if (rpX > rpU)
-                            colUp++;
-                          else
-                            {
-                              v -= retval.xdata (colXp) * data (colUp);
-                              colXp++;
-                              colUp++;
-                            }
-                        }
-                      while (rpX < j && rpU < j && colXp < cx && colUp < nz);
-
-                      // get A(m,m)
-                      if (typ == MatrixType::Upper)
-                        colUp = cidx (j+1) - 1;
-                      else
-                        colUp = cidx (j);
-                      double pivot = data (colUp);
-                      if (pivot == 0. || ridx (colUp) != j)
-                        {
-                          (*current_liboctave_error_handler)
-                            ("division by zero");
-                          goto inverse_singular;
-                        }
-
-                      if (v != 0.)
-                        {
-                          if (cx == nz2)
-                            {
-                              nz2 *= 2;
-                              retval.change_capacity (nz2);
-                            }
-
-                          retval.xridx (cx) = j;
-                          retval.xdata (cx) = v / pivot;
-                          cx++;
-                        }
-                    }
-
-                  // get A(m,m)
-                  octave_idx_type colUp;
-                  if (typ == MatrixType::Upper)
-                    colUp = cidx (i+1) - 1;
-                  else
-                    colUp = cidx (i);
-                  double pivot = data (colUp);
-                  if (pivot == 0. || ridx (colUp) != i)
-                    {
-                      (*current_liboctave_error_handler) ("division by zero");
-                      goto inverse_singular;
-                    }
-
-                  if (pivot != 1.0)
-                    for (octave_idx_type j = cx_colstart; j < cx; j++)
-                      retval.xdata (j) /= pivot;
                 }
-              retval.xcidx (nr) = cx;
-              retval.maybe_compress ();
-            }
+            }
+
+          // get A(m,m)
+          octave_idx_type colUp;
+          if (typ == MatrixType::Upper)
+            colUp = cidx (i+1) - 1;
           else
-            {
-              octave_idx_type nz = nnz ();
-              octave_idx_type cx = 0;
-              octave_idx_type nz2 = nz;
-              retval = SparseMatrix (nr, nc, nz2);
-
-              OCTAVE_LOCAL_BUFFER (double, work, nr);
-              OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr);
-
-              octave_idx_type *perm = mattyp.triangular_perm ();
-              if (typ == MatrixType::Permuted_Upper)
-                {
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    rperm[perm[i]] = i;
-                }
-              else
-                {
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    rperm[i] = perm[i];
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    perm[rperm[i]] = i;
-                }
-
-              for (octave_idx_type i = 0; i < nr; i++)
+            colUp = cidx (i);
+          double pivot = data (colUp);
+          if (pivot == 0. || ridx (colUp) != i)
+            (*current_liboctave_error_handler) ("division by zero");
+
+          if (pivot != 1.0)
+            for (octave_idx_type j = cx_colstart; j < cx; j++)
+              retval.xdata (j) /= pivot;
+        }
+      retval.xcidx (nr) = cx;
+      retval.maybe_compress ();
+    }
+  else
+    {
+      octave_idx_type nz = nnz ();
+      octave_idx_type cx = 0;
+      octave_idx_type nz2 = nz;
+      retval = SparseMatrix (nr, nc, nz2);
+
+      OCTAVE_LOCAL_BUFFER (double, work, nr);
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr);
+
+      octave_idx_type *perm = mattyp.triangular_perm ();
+      if (typ == MatrixType::Permuted_Upper)
+        {
+          for (octave_idx_type i = 0; i < nr; i++)
+            rperm[perm[i]] = i;
+        }
+      else
+        {
+          for (octave_idx_type i = 0; i < nr; i++)
+            rperm[i] = perm[i];
+          for (octave_idx_type i = 0; i < nr; i++)
+            perm[rperm[i]] = i;
+        }
+
+      for (octave_idx_type i = 0; i < nr; i++)
+        {
+          octave_quit ();
+          octave_idx_type iidx = rperm[i];
+
+          for (octave_idx_type j = 0; j < nr; j++)
+            work[j] = 0.;
+
+          // place the 1 in the identity position
+          work[iidx] = 1.0;
+
+          // iterate accross columns of input matrix
+          for (octave_idx_type j = iidx+1; j < nr; j++)
+            {
+              double v = 0.;
+              octave_idx_type jidx = perm[j];
+              // iterate to calculate sum
+              for (octave_idx_type k = cidx (jidx);
+                   k < cidx (jidx+1); k++)
                 {
                   octave_quit ();
-                  octave_idx_type iidx = rperm[i];
-
-                  for (octave_idx_type j = 0; j < nr; j++)
-                    work[j] = 0.;
-
-                  // place the 1 in the identity position
-                  work[iidx] = 1.0;
-
-                  // iterate accross columns of input matrix
-                  for (octave_idx_type j = iidx+1; j < nr; j++)
-                    {
-                      double v = 0.;
-                      octave_idx_type jidx = perm[j];
-                      // iterate to calculate sum
-                      for (octave_idx_type k = cidx (jidx);
-                           k < cidx (jidx+1); k++)
-                        {
-                          octave_quit ();
-                          v -= work[ridx (k)] * data (k);
-                        }
-
-                      // get A(m,m)
-                      double pivot;
-                      if (typ == MatrixType::Permuted_Upper)
-                        pivot = data (cidx (jidx+1) - 1);
-                      else
-                        pivot = data (cidx (jidx));
-                      if (pivot == 0.)
-                        {
-                          (*current_liboctave_error_handler)
-                            ("division by zero");
-                          goto inverse_singular;
-                        }
-
-                      work[j] = v / pivot;
-                    }
-
-                  // get A(m,m)
-                  octave_idx_type colUp;
-                  if (typ == MatrixType::Permuted_Upper)
-                    colUp = cidx (perm[iidx]+1) - 1;
-                  else
-                    colUp = cidx (perm[iidx]);
-
-                  double pivot = data (colUp);
-                  if (pivot == 0.)
-                    {
-                      (*current_liboctave_error_handler)
-                        ("division by zero");
-                      goto inverse_singular;
-                    }
-
-                  octave_idx_type new_cx = cx;
-                  for (octave_idx_type j = iidx; j < nr; j++)
-                    if (work[j] != 0.0)
-                      {
-                        new_cx++;
-                        if (pivot != 1.0)
-                          work[j] /= pivot;
-                      }
-
-                  if (cx < new_cx)
-                    {
-                      nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
-                      retval.change_capacity (nz2);
-                    }
-
-                  retval.xcidx (i) = cx;
-                  for (octave_idx_type j = iidx; j < nr; j++)
-                    if (work[j] != 0.)
-                      {
-                        retval.xridx (cx) = j;
-                        retval.xdata (cx++) = work[j];
-                      }
+                  v -= work[ridx (k)] * data (k);
                 }
 
-              retval.xcidx (nr) = cx;
-              retval.maybe_compress ();
-            }
-
-          if (calccond)
-            {
-              // Calculate the 1-norm of inverse matrix for rcond calculation
-              for (octave_idx_type j = 0; j < nr; j++)
-                {
-                  double atmp = 0.;
-                  for (octave_idx_type i = retval.cidx (j);
-                       i < retval.cidx (j+1); i++)
-                    atmp += fabs (retval.data (i));
-                  if (atmp > ainvnorm)
-                    ainvnorm = atmp;
-                }
-
-              rcond = 1. / ainvnorm / anorm;
-            }
+              // get A(m,m)
+              double pivot;
+              if (typ == MatrixType::Permuted_Upper)
+                pivot = data (cidx (jidx+1) - 1);
+              else
+                pivot = data (cidx (jidx));
+              if (pivot == 0.)
+                  (*current_liboctave_error_handler) ("division by zero");
+
+              work[j] = v / pivot;
+            }
+
+          // get A(m,m)
+          octave_idx_type colUp;
+          if (typ == MatrixType::Permuted_Upper)
+            colUp = cidx (perm[iidx]+1) - 1;
+          else
+            colUp = cidx (perm[iidx]);
+
+          double pivot = data (colUp);
+          if (pivot == 0.)
+            (*current_liboctave_error_handler) ("division by zero");
+
+          octave_idx_type new_cx = cx;
+          for (octave_idx_type j = iidx; j < nr; j++)
+            if (work[j] != 0.0)
+              {
+                new_cx++;
+                if (pivot != 1.0)
+                  work[j] /= pivot;
+              }
+
+          if (cx < new_cx)
+            {
+              nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
+              retval.change_capacity (nz2);
+            }
+
+          retval.xcidx (i) = cx;
+          for (octave_idx_type j = iidx; j < nr; j++)
+            if (work[j] != 0.)
+              {
+                retval.xridx (cx) = j;
+                retval.xdata (cx++) = work[j];
+              }
         }
-      else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      retval.xcidx (nr) = cx;
+      retval.maybe_compress ();
+    }
+
+  if (calccond)
+    {
+      // Calculate the 1-norm of inverse matrix for rcond calculation
+      for (octave_idx_type j = 0; j < nr; j++)
+        {
+          double atmp = 0.;
+          for (octave_idx_type i = retval.cidx (j);
+               i < retval.cidx (j+1); i++)
+            atmp += fabs (retval.data (i));
+          if (atmp > ainvnorm)
+            ainvnorm = atmp;
+        }
+
+      rcond = 1. / ainvnorm / anorm;
     }
 
   return retval;
@@ -1312,13 +1283,13 @@
 
       if (status < 0)
         {
-          (*current_liboctave_error_handler)
-            ("SparseMatrix::determinant symbolic factorization failed");
-
           UMFPACK_DNAME (report_status) (control, status);
           UMFPACK_DNAME (report_info) (control, info);
 
           UMFPACK_DNAME (free_symbolic) (&Symbolic);
+
+          (*current_liboctave_error_handler)
+            ("SparseMatrix::determinant symbolic factorization failed");
         }
       else
         {
@@ -1333,13 +1304,12 @@
 
           if (status < 0)
             {
-              (*current_liboctave_error_handler)
-                ("SparseMatrix::determinant numeric factorization failed");
-
               UMFPACK_DNAME (report_status) (control, status);
               UMFPACK_DNAME (report_info) (control, info);
 
               UMFPACK_DNAME (free_numeric) (&Numeric);
+              (*current_liboctave_error_handler)
+                ("SparseMatrix::determinant numeric factorization failed");
             }
           else
             {
@@ -1352,11 +1322,11 @@
 
               if (status < 0)
                 {
+                  UMFPACK_DNAME (report_status) (control, status);
+                  UMFPACK_DNAME (report_info) (control, info);
+
                   (*current_liboctave_error_handler)
                     ("SparseMatrix::determinant error calculating determinant");
-
-                  UMFPACK_DNAME (report_status) (control, status);
-                  UMFPACK_DNAME (report_info) (control, info);
                 }
               else
                 retval = DET (c10, e10, 10);
@@ -1391,7 +1361,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = Matrix (nc, b.cols (), 0.0);
   else
     {
@@ -1399,38 +1370,36 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
+      if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      retval.resize (nc, b.cols (), 0.);
+      if (typ == MatrixType::Diagonal)
+        for (octave_idx_type j = 0; j < b.cols (); j++)
+          for (octave_idx_type i = 0; i < nm; i++)
+            retval(i,j) = b(i,j) / data (i);
+      else
+        for (octave_idx_type j = 0; j < b.cols (); j++)
+          for (octave_idx_type k = 0; k < nc; k++)
+            for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+              retval(k,j) = b(ridx (i),j) / data (i);
+
+      if (calc_cond)
         {
-          retval.resize (nc, b.cols (), 0.);
-          if (typ == MatrixType::Diagonal)
-            for (octave_idx_type j = 0; j < b.cols (); j++)
-              for (octave_idx_type i = 0; i < nm; i++)
-                retval(i,j) = b(i,j) / data (i);
-          else
-            for (octave_idx_type j = 0; j < b.cols (); j++)
-              for (octave_idx_type k = 0; k < nc; k++)
-                for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
-                  retval(k,j) = b(ridx (i),j) / data (i);
-
-          if (calc_cond)
-            {
-              double dmax = 0.;
-              double dmin = octave_Inf;
-              for (octave_idx_type i = 0; i < nm; i++)
-                {
-                  double tmp = fabs (data (i));
-                  if (tmp > dmax)
-                    dmax = tmp;
-                  if (tmp < dmin)
-                    dmin = tmp;
-                }
-              rcond = dmin / dmax;
-            }
-          else
-            rcond = 1.;
+          double dmax = 0.;
+          double dmin = octave_Inf;
+          for (octave_idx_type i = 0; i < nm; i++)
+            {
+              double tmp = fabs (data (i));
+              if (tmp > dmax)
+                dmax = tmp;
+              if (tmp < dmin)
+                dmin = tmp;
+            }
+          rcond = dmin / dmax;
         }
       else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+        rcond = 1.;
     }
 
   return retval;
@@ -1451,7 +1420,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = SparseMatrix (nc, b.cols ());
   else
     {
@@ -1459,68 +1429,66 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
-        {
-          octave_idx_type b_nc = b.cols ();
-          octave_idx_type b_nz = b.nnz ();
-          retval = SparseMatrix (nc, b_nc, b_nz);
-
-          retval.xcidx (0) = 0;
-          octave_idx_type ii = 0;
-          if (typ == MatrixType::Diagonal)
-            for (octave_idx_type j = 0; j < b_nc; j++)
+      if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      octave_idx_type b_nc = b.cols ();
+      octave_idx_type b_nz = b.nnz ();
+      retval = SparseMatrix (nc, b_nc, b_nz);
+
+      retval.xcidx (0) = 0;
+      octave_idx_type ii = 0;
+      if (typ == MatrixType::Diagonal)
+        for (octave_idx_type j = 0; j < b_nc; j++)
+          {
+            for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
               {
-                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;
+                if (b.ridx (i) >= nm)
+                  break;
+                retval.xridx (ii) = b.ridx (i);
+                retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
               }
-          else
-            for (octave_idx_type j = 0; j < b_nc; j++)
-              {
-                for (octave_idx_type l = 0; l < nc; l++)
-                  for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
+            retval.xcidx (j+1) = ii;
+          }
+      else
+        for (octave_idx_type j = 0; j < b_nc; j++)
+          {
+            for (octave_idx_type l = 0; l < nc; l++)
+              for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
+                {
+                  bool found = false;
+                  octave_idx_type k;
+                  for (k = b.cidx (j); k < b.cidx (j+1); k++)
+                    if (ridx (i) == b.ridx (k))
+                      {
+                        found = true;
+                        break;
+                      }
+                  if (found)
                     {
-                      bool found = false;
-                      octave_idx_type k;
-                      for (k = b.cidx (j); k < b.cidx (j+1); k++)
-                        if (ridx (i) == b.ridx (k))
-                          {
-                            found = true;
-                            break;
-                          }
-                      if (found)
-                        {
-                          retval.xridx (ii) = l;
-                          retval.xdata (ii++) = b.data (k) / data (i);
-                        }
+                      retval.xridx (ii) = l;
+                      retval.xdata (ii++) = b.data (k) / data (i);
                     }
-                retval.xcidx (j+1) = ii;
-              }
-
-          if (calc_cond)
-            {
-              double dmax = 0.;
-              double dmin = octave_Inf;
-              for (octave_idx_type i = 0; i < nm; i++)
-                {
-                  double tmp = fabs (data (i));
-                  if (tmp > dmax)
-                    dmax = tmp;
-                  if (tmp < dmin)
-                    dmin = tmp;
                 }
-              rcond = dmin / dmax;
-            }
-          else
-            rcond = 1.;
+            retval.xcidx (j+1) = ii;
+          }
+
+      if (calc_cond)
+        {
+          double dmax = 0.;
+          double dmin = octave_Inf;
+          for (octave_idx_type i = 0; i < nm; i++)
+            {
+              double tmp = fabs (data (i));
+              if (tmp > dmax)
+                dmax = tmp;
+              if (tmp < dmin)
+                dmin = tmp;
+            }
+          rcond = dmin / dmax;
         }
       else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+        rcond = 1.;
     }
 
   return retval;
@@ -1541,7 +1509,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
   else
     {
@@ -1549,38 +1518,36 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
+      if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      retval.resize (nc, b.cols (), 0);
+      if (typ == MatrixType::Diagonal)
+        for (octave_idx_type j = 0; j < b.cols (); j++)
+          for (octave_idx_type i = 0; i < nm; i++)
+            retval(i,j) = b(i,j) / data (i);
+      else
+        for (octave_idx_type j = 0; j < b.cols (); j++)
+          for (octave_idx_type k = 0; k < nc; k++)
+            for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+              retval(k,j) = b(ridx (i),j) / data (i);
+
+      if (calc_cond)
         {
-          retval.resize (nc, b.cols (), 0);
-          if (typ == MatrixType::Diagonal)
-            for (octave_idx_type j = 0; j < b.cols (); j++)
-              for (octave_idx_type i = 0; i < nm; i++)
-                retval(i,j) = b(i,j) / data (i);
-          else
-            for (octave_idx_type j = 0; j < b.cols (); j++)
-              for (octave_idx_type k = 0; k < nc; k++)
-                for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
-                  retval(k,j) = b(ridx (i),j) / data (i);
-
-          if (calc_cond)
-            {
-              double dmax = 0.;
-              double dmin = octave_Inf;
-              for (octave_idx_type i = 0; i < nm; i++)
-                {
-                  double tmp = fabs (data (i));
-                  if (tmp > dmax)
-                    dmax = tmp;
-                  if (tmp < dmin)
-                    dmin = tmp;
-                }
-              rcond = dmin / dmax;
-            }
-          else
-            rcond = 1.;
+          double dmax = 0.;
+          double dmin = octave_Inf;
+          for (octave_idx_type i = 0; i < nm; i++)
+            {
+              double tmp = fabs (data (i));
+              if (tmp > dmax)
+                dmax = tmp;
+              if (tmp < dmin)
+                dmin = tmp;
+            }
+          rcond = dmin / dmax;
         }
       else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+        rcond = 1.;
     }
 
   return retval;
@@ -1601,7 +1568,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = SparseComplexMatrix (nc, b.cols ());
   else
     {
@@ -1609,68 +1577,66 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
-        {
-          octave_idx_type b_nc = b.cols ();
-          octave_idx_type b_nz = b.nnz ();
-          retval = SparseComplexMatrix (nc, b_nc, b_nz);
-
-          retval.xcidx (0) = 0;
-          octave_idx_type ii = 0;
-          if (typ == MatrixType::Diagonal)
-            for (octave_idx_type j = 0; j < b.cols (); j++)
+      if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      octave_idx_type b_nc = b.cols ();
+      octave_idx_type b_nz = b.nnz ();
+      retval = SparseComplexMatrix (nc, b_nc, b_nz);
+
+      retval.xcidx (0) = 0;
+      octave_idx_type ii = 0;
+      if (typ == MatrixType::Diagonal)
+        for (octave_idx_type j = 0; j < b.cols (); j++)
+          {
+            for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
               {
-                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;
+                if (b.ridx (i) >= nm)
+                  break;
+                retval.xridx (ii) = b.ridx (i);
+                retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
               }
-          else
-            for (octave_idx_type j = 0; j < b.cols (); j++)
-              {
-                for (octave_idx_type l = 0; l < nc; l++)
-                  for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
+            retval.xcidx (j+1) = ii;
+          }
+      else
+        for (octave_idx_type j = 0; j < b.cols (); j++)
+          {
+            for (octave_idx_type l = 0; l < nc; l++)
+              for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
+                {
+                  bool found = false;
+                  octave_idx_type k;
+                  for (k = b.cidx (j); k < b.cidx (j+1); k++)
+                    if (ridx (i) == b.ridx (k))
+                      {
+                        found = true;
+                        break;
+                      }
+                  if (found)
                     {
-                      bool found = false;
-                      octave_idx_type k;
-                      for (k = b.cidx (j); k < b.cidx (j+1); k++)
-                        if (ridx (i) == b.ridx (k))
-                          {
-                            found = true;
-                            break;
-                          }
-                      if (found)
-                        {
-                          retval.xridx (ii) = l;
-                          retval.xdata (ii++) = b.data (k) / data (i);
-                        }
+                      retval.xridx (ii) = l;
+                      retval.xdata (ii++) = b.data (k) / data (i);
                     }
-                retval.xcidx (j+1) = ii;
-              }
-
-          if (calc_cond)
-            {
-              double dmax = 0.;
-              double dmin = octave_Inf;
-              for (octave_idx_type i = 0; i < nm; i++)
-                {
-                  double tmp = fabs (data (i));
-                  if (tmp > dmax)
-                    dmax = tmp;
-                  if (tmp < dmin)
-                    dmin = tmp;
                 }
-              rcond = dmin / dmax;
-            }
-          else
-            rcond = 1.;
+            retval.xcidx (j+1) = ii;
+          }
+
+      if (calc_cond)
+        {
+          double dmax = 0.;
+          double dmin = octave_Inf;
+          for (octave_idx_type i = 0; i < nm; i++)
+            {
+              double tmp = fabs (data (i));
+              if (tmp > dmax)
+                dmax = tmp;
+              if (tmp < dmin)
+                dmin = tmp;
+            }
+          rcond = dmin / dmax;
         }
       else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+        rcond = 1.;
     }
 
   return retval;
@@ -1692,7 +1658,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = Matrix (nc, b.cols (), 0.0);
   else
     {
@@ -1700,128 +1667,157 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper)
+      if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      double anorm = 0.;
+      double ainvnorm = 0.;
+      octave_idx_type b_nc = b.cols ();
+      rcond = 1.;
+
+      if (calc_cond)
+        {
+          // Calculate the 1-norm of matrix for rcond calculation
+          for (octave_idx_type j = 0; j < nc; j++)
+            {
+              double atmp = 0.;
+              for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                atmp += fabs (data (i));
+              if (atmp > anorm)
+                anorm = atmp;
+            }
+        }
+
+      if (typ == MatrixType::Permuted_Upper)
         {
-          double anorm = 0.;
-          double ainvnorm = 0.;
-          octave_idx_type b_nc = b.cols ();
-          rcond = 1.;
+          retval.resize (nc, b_nc);
+          octave_idx_type *perm = mattype.triangular_perm ();
+          OCTAVE_LOCAL_BUFFER (double, work, nm);
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nr; i++)
+                work[i] = b(i,j);
+              for (octave_idx_type i = nr; i < nc; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type k = nc-1; k >= 0; k--)
+                {
+                  octave_idx_type kidx = perm[k];
+
+                  if (work[k] != 0.)
+                    {
+                      if (ridx (cidx (kidx+1)-1) != k
+                          || data (cidx (kidx+1)-1) == 0.)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      double tmp = work[k] / data (cidx (kidx+1)-1);
+                      work[k] = tmp;
+                      for (octave_idx_type i = cidx (kidx);
+                           i < cidx (kidx+1)-1; i++)
+                        {
+                          octave_idx_type iidx = ridx (i);
+                          work[iidx] = work[iidx] - tmp * data (i);
+                        }
+                    }
+                }
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                retval.xelem (perm[i], j) = work[i];
+            }
 
           if (calc_cond)
             {
-              // Calculate the 1-norm of matrix for rcond calculation
-              for (octave_idx_type j = 0; j < nc; j++)
+              // Calculation of 1-norm of inv(*this)
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  double atmp = 0.;
-                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
-                    atmp += fabs (data (i));
-                  if (atmp > anorm)
-                    anorm = atmp;
-                }
-            }
-
-          if (typ == MatrixType::Permuted_Upper)
-            {
-              retval.resize (nc, b_nc);
-              octave_idx_type *perm = mattype.triangular_perm ();
-              OCTAVE_LOCAL_BUFFER (double, work, nm);
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
-                {
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    work[i] = b(i,j);
-                  for (octave_idx_type i = nr; i < nc; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k >= 0; k--)
                     {
-                      octave_idx_type kidx = perm[k];
+                      octave_idx_type iidx = perm[k];
 
                       if (work[k] != 0.)
                         {
-                          if (ridx (cidx (kidx+1)-1) != k
-                              || data (cidx (kidx+1)-1) == 0.)
+                          double tmp = work[k] / data (cidx (iidx+1)-1);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (iidx);
+                               i < cidx (iidx+1)-1; i++)
                             {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
-                          double tmp = work[k] / data (cidx (kidx+1)-1);
-                          work[k] = tmp;
-                          for (octave_idx_type i = cidx (kidx);
-                               i < cidx (kidx+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
+                              octave_idx_type idx2 = ridx (i);
+                              work[idx2] = work[idx2] - tmp * data (i);
                             }
                         }
                     }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    retval.xelem (perm[i], j) = work[i];
+                  double atmp = 0;
+                  for (octave_idx_type i = 0; i < j+1; i++)
+                    {
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
+                    }
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-
-              if (calc_cond)
+              rcond = 1. / ainvnorm / anorm;
+            }
+        }
+      else
+        {
+          OCTAVE_LOCAL_BUFFER (double, work, nm);
+          retval.resize (nc, b_nc);
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nr; i++)
+                work[i] = b(i,j);
+              for (octave_idx_type i = nr; i < nc; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type k = nc-1; k >= 0; k--)
                 {
-                  // Calculation of 1-norm of inv(*this)
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  if (work[k] != 0.)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k >= 0; k--)
+                      if (ridx (cidx (k+1)-1) != k
+                          || data (cidx (k+1)-1) == 0.)
                         {
-                          octave_idx_type iidx = perm[k];
-
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (iidx+1)-1);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (iidx);
-                                   i < cidx (iidx+1)-1; i++)
-                                {
-                                  octave_idx_type idx2 = ridx (i);
-                                  work[idx2] = work[idx2] - tmp * data (i);
-                                }
-                            }
+                          err = -2;
+                          goto triangular_error;
                         }
-                      double atmp = 0;
-                      for (octave_idx_type i = 0; i < j+1; i++)
+
+                      double tmp = work[k] / data (cidx (k+1)-1);
+                      work[k] = tmp;
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
                         {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
+                          octave_idx_type iidx = ridx (i);
+                          work[iidx] = work[iidx] - tmp * data (i);
                         }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
                     }
-                  rcond = 1. / ainvnorm / anorm;
                 }
-            }
-          else
-            {
-              OCTAVE_LOCAL_BUFFER (double, work, nm);
-              retval.resize (nc, b_nc);
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                retval.xelem (i, j) = work[i];
+            }
+
+          if (calc_cond)
+            {
+              // Calculation of 1-norm of inv(*this)
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    work[i] = b(i,j);
-                  for (octave_idx_type i = nr; i < nc; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k >= 0; k--)
                     {
                       if (work[k] != 0.)
                         {
-                          if (ridx (cidx (k+1)-1) != k
-                              || data (cidx (k+1)-1) == 0.)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
                           double tmp = work[k] / data (cidx (k+1)-1);
                           work[k] = tmp;
                           for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
@@ -1831,76 +1827,45 @@
                             }
                         }
                     }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    retval.xelem (i, j) = work[i];
-                }
-
-              if (calc_cond)
-                {
-                  // Calculation of 1-norm of inv(*this)
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  double atmp = 0;
+                  for (octave_idx_type i = 0; i < j+1; i++)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k >= 0; k--)
-                        {
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (k+1)-1);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
-                        }
-                      double atmp = 0;
-                      for (octave_idx_type i = 0; i < j+1; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
                     }
-                  rcond = 1. / ainvnorm / anorm;
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-            }
-
-        triangular_error:
-          if (err != 0)
-            {
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
-            }
-
-          volatile double rcond_plus_one = rcond + 1.0;
-
-          if (rcond_plus_one == 1.0 || xisnan (rcond))
-            {
-              err = -2;
-
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
+              rcond = 1. / ainvnorm / anorm;
             }
         }
-      else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+    triangular_error:
+      if (err != 0)
+        {
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
+
+      volatile double rcond_plus_one = rcond + 1.0;
+
+      if (rcond_plus_one == 1.0 || xisnan (rcond))
+        {
+          err = -2;
+
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
     }
 
   return retval;
@@ -1922,7 +1887,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = SparseMatrix (nc, b.cols ());
   else
     {
@@ -1930,260 +1896,258 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper)
+      if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      double anorm = 0.;
+      double ainvnorm = 0.;
+      rcond = 1.;
+
+      if (calc_cond)
+        {
+          // Calculate the 1-norm of matrix for rcond calculation
+          for (octave_idx_type j = 0; j < nc; j++)
+            {
+              double atmp = 0.;
+              for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                atmp += fabs (data (i));
+              if (atmp > anorm)
+                anorm = atmp;
+            }
+        }
+
+      octave_idx_type b_nc = b.cols ();
+      octave_idx_type b_nz = b.nnz ();
+      retval = SparseMatrix (nc, b_nc, b_nz);
+      retval.xcidx (0) = 0;
+      octave_idx_type ii = 0;
+      octave_idx_type x_nz = b_nz;
+
+      if (typ == MatrixType::Permuted_Upper)
         {
-          double anorm = 0.;
-          double ainvnorm = 0.;
-          rcond = 1.;
+          octave_idx_type *perm = mattype.triangular_perm ();
+          OCTAVE_LOCAL_BUFFER (double, work, nm);
+
+          OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
+          for (octave_idx_type i = 0; i < nc; i++)
+            rperm[perm[i]] = i;
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+              for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                work[b.ridx (i)] = b.data (i);
+
+              for (octave_idx_type k = nc-1; k >= 0; k--)
+                {
+                  octave_idx_type kidx = perm[k];
+
+                  if (work[k] != 0.)
+                    {
+                      if (ridx (cidx (kidx+1)-1) != k
+                          || data (cidx (kidx+1)-1) == 0.)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      double tmp = work[k] / data (cidx (kidx+1)-1);
+                      work[k] = tmp;
+                      for (octave_idx_type i = cidx (kidx);
+                           i < cidx (kidx+1)-1; i++)
+                        {
+                          octave_idx_type iidx = ridx (i);
+                          work[iidx] = work[iidx] - tmp * data (i);
+                        }
+                    }
+                }
+
+              // Count nonzeros in work vector and adjust space in
+              // retval if needed
+              octave_idx_type new_nnz = 0;
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (work[i] != 0.)
+                  new_nnz++;
+
+              if (ii + new_nnz > x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                  retval.change_capacity (sz);
+                  x_nz = sz;
+                }
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (work[rperm[i]] != 0.)
+                  {
+                    retval.xridx (ii) = i;
+                    retval.xdata (ii++) = work[rperm[i]];
+                  }
+              retval.xcidx (j+1) = ii;
+            }
+
+          retval.maybe_compress ();
 
           if (calc_cond)
             {
-              // Calculate the 1-norm of matrix for rcond calculation
-              for (octave_idx_type j = 0; j < nc; j++)
+              // Calculation of 1-norm of inv(*this)
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  double atmp = 0.;
-                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
-                    atmp += fabs (data (i));
-                  if (atmp > anorm)
-                    anorm = atmp;
-                }
-            }
-
-          octave_idx_type b_nc = b.cols ();
-          octave_idx_type b_nz = b.nnz ();
-          retval = SparseMatrix (nc, b_nc, b_nz);
-          retval.xcidx (0) = 0;
-          octave_idx_type ii = 0;
-          octave_idx_type x_nz = b_nz;
-
-          if (typ == MatrixType::Permuted_Upper)
-            {
-              octave_idx_type *perm = mattype.triangular_perm ();
-              OCTAVE_LOCAL_BUFFER (double, work, nm);
-
-              OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
-              for (octave_idx_type i = 0; i < nc; i++)
-                rperm[perm[i]] = i;
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
-                {
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
-                    work[b.ridx (i)] = b.data (i);
-
-                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k >= 0; k--)
                     {
-                      octave_idx_type kidx = perm[k];
+                      octave_idx_type iidx = perm[k];
 
                       if (work[k] != 0.)
                         {
-                          if (ridx (cidx (kidx+1)-1) != k
-                              || data (cidx (kidx+1)-1) == 0.)
+                          double tmp = work[k] / data (cidx (iidx+1)-1);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (iidx);
+                               i < cidx (iidx+1)-1; i++)
                             {
-                              err = -2;
-                              goto triangular_error;
+                              octave_idx_type idx2 = ridx (i);
+                              work[idx2] = work[idx2] - tmp * data (i);
                             }
-
-                          double tmp = work[k] / data (cidx (kidx+1)-1);
+                        }
+                    }
+                  double atmp = 0;
+                  for (octave_idx_type i = 0; i < j+1; i++)
+                    {
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
+                    }
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
+                }
+              rcond = 1. / ainvnorm / anorm;
+            }
+        }
+      else
+        {
+          OCTAVE_LOCAL_BUFFER (double, work, nm);
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+              for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                work[b.ridx (i)] = b.data (i);
+
+              for (octave_idx_type k = nc-1; k >= 0; k--)
+                {
+                  if (work[k] != 0.)
+                    {
+                      if (ridx (cidx (k+1)-1) != k
+                          || data (cidx (k+1)-1) == 0.)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      double tmp = work[k] / data (cidx (k+1)-1);
+                      work[k] = tmp;
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                        {
+                          octave_idx_type iidx = ridx (i);
+                          work[iidx] = work[iidx] - tmp * data (i);
+                        }
+                    }
+                }
+
+              // Count nonzeros in work vector and adjust space in
+              // retval if needed
+              octave_idx_type new_nnz = 0;
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (work[i] != 0.)
+                  new_nnz++;
+
+              if (ii + new_nnz > x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                  retval.change_capacity (sz);
+                  x_nz = sz;
+                }
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (work[i] != 0.)
+                  {
+                    retval.xridx (ii) = i;
+                    retval.xdata (ii++) = work[i];
+                  }
+              retval.xcidx (j+1) = ii;
+            }
+
+          retval.maybe_compress ();
+
+          if (calc_cond)
+            {
+              // Calculation of 1-norm of inv(*this)
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
+                {
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k >= 0; k--)
+                    {
+                      if (work[k] != 0.)
+                        {
+                          double tmp = work[k] / data (cidx (k+1)-1);
                           work[k] = tmp;
-                          for (octave_idx_type i = cidx (kidx);
-                               i < cidx (kidx+1)-1; i++)
+                          for (octave_idx_type i = cidx (k);
+                               i < cidx (k+1)-1; i++)
                             {
                               octave_idx_type iidx = ridx (i);
                               work[iidx] = work[iidx] - tmp * data (i);
                             }
                         }
                     }
-
-                  // Count nonzeros in work vector and adjust space in
-                  // retval if needed
-                  octave_idx_type new_nnz = 0;
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (work[i] != 0.)
-                      new_nnz++;
-
-                  if (ii + new_nnz > x_nz)
+                  double atmp = 0;
+                  for (octave_idx_type i = 0; i < j+1; i++)
                     {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
-                      retval.change_capacity (sz);
-                      x_nz = sz;
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
                     }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (work[rperm[i]] != 0.)
-                      {
-                        retval.xridx (ii) = i;
-                        retval.xdata (ii++) = work[rperm[i]];
-                      }
-                  retval.xcidx (j+1) = ii;
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-
-              retval.maybe_compress ();
-
-              if (calc_cond)
-                {
-                  // Calculation of 1-norm of inv(*this)
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
-                    {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k >= 0; k--)
-                        {
-                          octave_idx_type iidx = perm[k];
-
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (iidx+1)-1);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (iidx);
-                                   i < cidx (iidx+1)-1; i++)
-                                {
-                                  octave_idx_type idx2 = ridx (i);
-                                  work[idx2] = work[idx2] - tmp * data (i);
-                                }
-                            }
-                        }
-                      double atmp = 0;
-                      for (octave_idx_type i = 0; i < j+1; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
-                    }
-                  rcond = 1. / ainvnorm / anorm;
-                }
+              rcond = 1. / ainvnorm / anorm;
+            }
+        }
+
+    triangular_error:
+      if (err != 0)
+        {
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
             }
           else
-            {
-              OCTAVE_LOCAL_BUFFER (double, work, nm);
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
-                {
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
-                    work[b.ridx (i)] = b.data (i);
-
-                  for (octave_idx_type k = nc-1; k >= 0; k--)
-                    {
-                      if (work[k] != 0.)
-                        {
-                          if (ridx (cidx (k+1)-1) != k
-                              || data (cidx (k+1)-1) == 0.)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
-                          double tmp = work[k] / data (cidx (k+1)-1);
-                          work[k] = tmp;
-                          for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              work[iidx] = work[iidx] - tmp * data (i);
-                            }
-                        }
-                    }
-
-                  // Count nonzeros in work vector and adjust space in
-                  // retval if needed
-                  octave_idx_type new_nnz = 0;
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (work[i] != 0.)
-                      new_nnz++;
-
-                  if (ii + new_nnz > x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
-                      retval.change_capacity (sz);
-                      x_nz = sz;
-                    }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (work[i] != 0.)
-                      {
-                        retval.xridx (ii) = i;
-                        retval.xdata (ii++) = work[i];
-                      }
-                  retval.xcidx (j+1) = ii;
-                }
-
-              retval.maybe_compress ();
-
-              if (calc_cond)
-                {
-                  // Calculation of 1-norm of inv(*this)
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
-                    {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k >= 0; k--)
-                        {
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (k+1)-1);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1)-1; i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
-                        }
-                      double atmp = 0;
-                      for (octave_idx_type i = 0; i < j+1; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
-                    }
-                  rcond = 1. / ainvnorm / anorm;
-                }
-            }
-
-        triangular_error:
-          if (err != 0)
-            {
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
-            }
-
-          volatile double rcond_plus_one = rcond + 1.0;
-
-          if (rcond_plus_one == 1.0 || xisnan (rcond))
-            {
-              err = -2;
-
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
-            }
+            warn_singular_matrix (rcond);
         }
-      else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      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
+            warn_singular_matrix (rcond);
+        }
     }
   return retval;
 }
@@ -2204,7 +2168,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
   else
     {
@@ -2212,210 +2177,208 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper)
+      if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      double anorm = 0.;
+      double ainvnorm = 0.;
+      octave_idx_type b_nc = b.cols ();
+      rcond = 1.;
+
+      if (calc_cond)
+        {
+          // Calculate the 1-norm of matrix for rcond calculation
+          for (octave_idx_type j = 0; j < nc; j++)
+            {
+              double atmp = 0.;
+              for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                atmp += fabs (data (i));
+              if (atmp > anorm)
+                anorm = atmp;
+            }
+        }
+
+      if (typ == MatrixType::Permuted_Upper)
         {
-          double anorm = 0.;
-          double ainvnorm = 0.;
-          octave_idx_type b_nc = b.cols ();
-          rcond = 1.;
+          retval.resize (nc, b_nc);
+          octave_idx_type *perm = mattype.triangular_perm ();
+          OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nr; i++)
+                cwork[i] = b(i,j);
+              for (octave_idx_type i = nr; i < nc; i++)
+                cwork[i] = 0.;
+
+              for (octave_idx_type k = nc-1; k >= 0; k--)
+                {
+                  octave_idx_type kidx = perm[k];
+
+                  if (cwork[k] != 0.)
+                    {
+                      if (ridx (cidx (kidx+1)-1) != k
+                          || data (cidx (kidx+1)-1) == 0.)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
+                      cwork[k] = tmp;
+                      for (octave_idx_type i = cidx (kidx);
+                           i < cidx (kidx+1)-1; i++)
+                        {
+                          octave_idx_type iidx = ridx (i);
+                          cwork[iidx] = cwork[iidx] - tmp * data (i);
+                        }
+                    }
+                }
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                retval.xelem (perm[i], j) = cwork[i];
+            }
 
           if (calc_cond)
             {
-              // Calculate the 1-norm of matrix for rcond calculation
-              for (octave_idx_type j = 0; j < nc; j++)
-                {
-                  double atmp = 0.;
-                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
-                    atmp += fabs (data (i));
-                  if (atmp > anorm)
-                    anorm = atmp;
-                }
-            }
-
-          if (typ == MatrixType::Permuted_Upper)
-            {
-              retval.resize (nc, b_nc);
-              octave_idx_type *perm = mattype.triangular_perm ();
-              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
+              // Calculation of 1-norm of inv(*this)
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    cwork[i] = b(i,j);
-                  for (octave_idx_type i = nr; i < nc; i++)
-                    cwork[i] = 0.;
-
-                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k >= 0; k--)
                     {
-                      octave_idx_type kidx = perm[k];
-
-                      if (cwork[k] != 0.)
+                      octave_idx_type iidx = perm[k];
+
+                      if (work[k] != 0.)
                         {
-                          if (ridx (cidx (kidx+1)-1) != k
-                              || data (cidx (kidx+1)-1) == 0.)
+                          double tmp = work[k] / data (cidx (iidx+1)-1);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (iidx);
+                               i < cidx (iidx+1)-1; i++)
                             {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
-                          Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
-                          cwork[k] = tmp;
-                          for (octave_idx_type i = cidx (kidx);
-                               i < cidx (kidx+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                              octave_idx_type idx2 = ridx (i);
+                              work[idx2] = work[idx2] - tmp * data (i);
                             }
                         }
                     }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    retval.xelem (perm[i], j) = cwork[i];
-                }
-
-              if (calc_cond)
-                {
-                  // Calculation of 1-norm of inv(*this)
-                  OCTAVE_LOCAL_BUFFER (double, work, nm);
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  double atmp = 0;
+                  for (octave_idx_type i = 0; i < j+1; i++)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k >= 0; k--)
-                        {
-                          octave_idx_type iidx = perm[k];
-
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (iidx+1)-1);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (iidx);
-                                   i < cidx (iidx+1)-1; i++)
-                                {
-                                  octave_idx_type idx2 = ridx (i);
-                                  work[idx2] = work[idx2] - tmp * data (i);
-                                }
-                            }
-                        }
-                      double atmp = 0;
-                      for (octave_idx_type i = 0; i < j+1; i++)
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
+                    }
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
+                }
+              rcond = 1. / ainvnorm / anorm;
+            }
+        }
+      else
+        {
+          OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+          retval.resize (nc, b_nc);
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nr; i++)
+                cwork[i] = b(i,j);
+              for (octave_idx_type i = nr; i < nc; i++)
+                cwork[i] = 0.;
+
+              for (octave_idx_type k = nc-1; k >= 0; k--)
+                {
+                  if (cwork[k] != 0.)
+                    {
+                      if (ridx (cidx (k+1)-1) != k
+                          || data (cidx (k+1)-1) == 0.)
                         {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
+                          err = -2;
+                          goto triangular_error;
                         }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
+
+                      Complex tmp = cwork[k] / data (cidx (k+1)-1);
+                      cwork[k] = tmp;
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                        {
+                          octave_idx_type iidx = ridx (i);
+                          cwork[iidx] = cwork[iidx] - tmp  * data (i);
+                        }
                     }
-                  rcond = 1. / ainvnorm / anorm;
                 }
-            }
-          else
-            {
-              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
-              retval.resize (nc, b_nc);
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                retval.xelem (i, j) = cwork[i];
+            }
+
+          if (calc_cond)
+            {
+              // Calculation of 1-norm of inv(*this)
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    cwork[i] = b(i,j);
-                  for (octave_idx_type i = nr; i < nc; i++)
-                    cwork[i] = 0.;
-
-                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k >= 0; k--)
                     {
-                      if (cwork[k] != 0.)
+                      if (work[k] != 0.)
                         {
-                          if (ridx (cidx (k+1)-1) != k
-                              || data (cidx (k+1)-1) == 0.)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
-                          Complex tmp = cwork[k] / data (cidx (k+1)-1);
-                          cwork[k] = tmp;
-                          for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                          double tmp = work[k] / data (cidx (k+1)-1);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k);
+                               i < cidx (k+1)-1; i++)
                             {
                               octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp  * data (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
                             }
                         }
                     }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    retval.xelem (i, j) = cwork[i];
-                }
-
-              if (calc_cond)
-                {
-                  // Calculation of 1-norm of inv(*this)
-                  OCTAVE_LOCAL_BUFFER (double, work, nm);
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  double atmp = 0;
+                  for (octave_idx_type i = 0; i < j+1; i++)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k >= 0; k--)
-                        {
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (k+1)-1);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1)-1; i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
-                        }
-                      double atmp = 0;
-                      for (octave_idx_type i = 0; i < j+1; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
                     }
-                  rcond = 1. / ainvnorm / anorm;
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-            }
-
-        triangular_error:
-          if (err != 0)
-            {
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
-            }
-
-          volatile double rcond_plus_one = rcond + 1.0;
-
-          if (rcond_plus_one == 1.0 || xisnan (rcond))
-            {
-              err = -2;
-
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
+              rcond = 1. / ainvnorm / anorm;
             }
         }
-      else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+    triangular_error:
+      if (err != 0)
+        {
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
+
+      volatile double rcond_plus_one = rcond + 1.0;
+
+      if (rcond_plus_one == 1.0 || xisnan (rcond))
+        {
+          err = -2;
+
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
     }
 
   return retval;
@@ -2437,7 +2400,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = SparseComplexMatrix (nc, b.cols ());
   else
     {
@@ -2445,262 +2409,260 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper)
+      if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      double anorm = 0.;
+      double ainvnorm = 0.;
+      rcond = 1.;
+
+      if (calc_cond)
+        {
+          // Calculate the 1-norm of matrix for rcond calculation
+          for (octave_idx_type j = 0; j < nc; j++)
+            {
+              double atmp = 0.;
+              for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                atmp += fabs (data (i));
+              if (atmp > anorm)
+                anorm = atmp;
+            }
+        }
+
+      octave_idx_type b_nc = b.cols ();
+      octave_idx_type b_nz = b.nnz ();
+      retval = SparseComplexMatrix (nc, b_nc, b_nz);
+      retval.xcidx (0) = 0;
+      octave_idx_type ii = 0;
+      octave_idx_type x_nz = b_nz;
+
+      if (typ == MatrixType::Permuted_Upper)
         {
-          double anorm = 0.;
-          double ainvnorm = 0.;
-          rcond = 1.;
+          octave_idx_type *perm = mattype.triangular_perm ();
+          OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+
+          OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
+          for (octave_idx_type i = 0; i < nc; i++)
+            rperm[perm[i]] = i;
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nm; i++)
+                cwork[i] = 0.;
+              for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                cwork[b.ridx (i)] = b.data (i);
+
+              for (octave_idx_type k = nc-1; k >= 0; k--)
+                {
+                  octave_idx_type kidx = perm[k];
+
+                  if (cwork[k] != 0.)
+                    {
+                      if (ridx (cidx (kidx+1)-1) != k
+                          || data (cidx (kidx+1)-1) == 0.)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
+                      cwork[k] = tmp;
+                      for (octave_idx_type i = cidx (kidx);
+                           i < cidx (kidx+1)-1; i++)
+                        {
+                          octave_idx_type iidx = ridx (i);
+                          cwork[iidx] = cwork[iidx] - tmp * data (i);
+                        }
+                    }
+                }
+
+              // Count nonzeros in work vector and adjust space in
+              // retval if needed
+              octave_idx_type new_nnz = 0;
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (cwork[i] != 0.)
+                  new_nnz++;
+
+              if (ii + new_nnz > x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                  retval.change_capacity (sz);
+                  x_nz = sz;
+                }
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (cwork[rperm[i]] != 0.)
+                  {
+                    retval.xridx (ii) = i;
+                    retval.xdata (ii++) = cwork[rperm[i]];
+                  }
+              retval.xcidx (j+1) = ii;
+            }
+
+          retval.maybe_compress ();
 
           if (calc_cond)
             {
-              // Calculate the 1-norm of matrix for rcond calculation
-              for (octave_idx_type j = 0; j < nc; j++)
+              // Calculation of 1-norm of inv(*this)
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  double atmp = 0.;
-                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
-                    atmp += fabs (data (i));
-                  if (atmp > anorm)
-                    anorm = atmp;
-                }
-            }
-
-          octave_idx_type b_nc = b.cols ();
-          octave_idx_type b_nz = b.nnz ();
-          retval = SparseComplexMatrix (nc, b_nc, b_nz);
-          retval.xcidx (0) = 0;
-          octave_idx_type ii = 0;
-          octave_idx_type x_nz = b_nz;
-
-          if (typ == MatrixType::Permuted_Upper)
-            {
-              octave_idx_type *perm = mattype.triangular_perm ();
-              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
-
-              OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
-              for (octave_idx_type i = 0; i < nc; i++)
-                rperm[perm[i]] = i;
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
-                {
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    cwork[i] = 0.;
-                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
-                    cwork[b.ridx (i)] = b.data (i);
-
-                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k >= 0; k--)
                     {
-                      octave_idx_type kidx = perm[k];
-
-                      if (cwork[k] != 0.)
+                      octave_idx_type iidx = perm[k];
+
+                      if (work[k] != 0.)
                         {
-                          if (ridx (cidx (kidx+1)-1) != k
-                              || data (cidx (kidx+1)-1) == 0.)
+                          double tmp = work[k] / data (cidx (iidx+1)-1);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (iidx);
+                               i < cidx (iidx+1)-1; i++)
                             {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
-                          Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
-                          cwork[k] = tmp;
-                          for (octave_idx_type i = cidx (kidx);
-                               i < cidx (kidx+1)-1; i++)
-                            {
-                              octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                              octave_idx_type idx2 = ridx (i);
+                              work[idx2] = work[idx2] - tmp * data (i);
                             }
                         }
                     }
-
-                  // Count nonzeros in work vector and adjust space in
-                  // retval if needed
-                  octave_idx_type new_nnz = 0;
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (cwork[i] != 0.)
-                      new_nnz++;
-
-                  if (ii + new_nnz > x_nz)
+                  double atmp = 0;
+                  for (octave_idx_type i = 0; i < j+1; i++)
                     {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
-                      retval.change_capacity (sz);
-                      x_nz = sz;
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
                     }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (cwork[rperm[i]] != 0.)
-                      {
-                        retval.xridx (ii) = i;
-                        retval.xdata (ii++) = cwork[rperm[i]];
-                      }
-                  retval.xcidx (j+1) = ii;
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-
-              retval.maybe_compress ();
-
-              if (calc_cond)
+              rcond = 1. / ainvnorm / anorm;
+            }
+        }
+      else
+        {
+          OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nm; i++)
+                cwork[i] = 0.;
+              for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                cwork[b.ridx (i)] = b.data (i);
+
+              for (octave_idx_type k = nc-1; k >= 0; k--)
                 {
-                  // Calculation of 1-norm of inv(*this)
-                  OCTAVE_LOCAL_BUFFER (double, work, nm);
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  if (cwork[k] != 0.)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k >= 0; k--)
+                      if (ridx (cidx (k+1)-1) != k
+                          || data (cidx (k+1)-1) == 0.)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      Complex tmp = cwork[k] / data (cidx (k+1)-1);
+                      cwork[k] = tmp;
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
                         {
-                          octave_idx_type iidx = perm[k];
-
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (iidx+1)-1);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (iidx);
-                                   i < cidx (iidx+1)-1; i++)
-                                {
-                                  octave_idx_type idx2 = ridx (i);
-                                  work[idx2] = work[idx2] - tmp * data (i);
-                                }
-                            }
+                          octave_idx_type iidx = ridx (i);
+                          cwork[iidx] = cwork[iidx] - tmp * data (i);
                         }
-                      double atmp = 0;
-                      for (octave_idx_type i = 0; i < j+1; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
                     }
-                  rcond = 1. / ainvnorm / anorm;
+                }
+
+              // Count nonzeros in work vector and adjust space in
+              // retval if needed
+              octave_idx_type new_nnz = 0;
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (cwork[i] != 0.)
+                  new_nnz++;
+
+              if (ii + new_nnz > x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                  retval.change_capacity (sz);
+                  x_nz = sz;
                 }
-            }
-          else
-            {
-              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (cwork[i] != 0.)
+                  {
+                    retval.xridx (ii) = i;
+                    retval.xdata (ii++) = cwork[i];
+                  }
+              retval.xcidx (j+1) = ii;
+            }
+
+          retval.maybe_compress ();
+
+          if (calc_cond)
+            {
+              // Calculation of 1-norm of inv(*this)
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    cwork[i] = 0.;
-                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
-                    cwork[b.ridx (i)] = b.data (i);
-
-                  for (octave_idx_type k = nc-1; k >= 0; k--)
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k >= 0; k--)
                     {
-                      if (cwork[k] != 0.)
+                      if (work[k] != 0.)
                         {
-                          if (ridx (cidx (k+1)-1) != k
-                              || data (cidx (k+1)-1) == 0.)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
-                          Complex tmp = cwork[k] / data (cidx (k+1)-1);
-                          cwork[k] = tmp;
-                          for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
+                          double tmp = work[k] / data (cidx (k+1)-1);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k);
+                               i < cidx (k+1)-1; i++)
                             {
                               octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
                             }
                         }
                     }
-
-                  // Count nonzeros in work vector and adjust space in
-                  // retval if needed
-                  octave_idx_type new_nnz = 0;
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (cwork[i] != 0.)
-                      new_nnz++;
-
-                  if (ii + new_nnz > x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
-                      retval.change_capacity (sz);
-                      x_nz = sz;
-                    }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (cwork[i] != 0.)
-                      {
-                        retval.xridx (ii) = i;
-                        retval.xdata (ii++) = cwork[i];
-                      }
-                  retval.xcidx (j+1) = ii;
-                }
-
-              retval.maybe_compress ();
-
-              if (calc_cond)
-                {
-                  // Calculation of 1-norm of inv(*this)
-                  OCTAVE_LOCAL_BUFFER (double, work, nm);
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  double atmp = 0;
+                  for (octave_idx_type i = 0; i < j+1; i++)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k >= 0; k--)
-                        {
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (k+1)-1);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1)-1; i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
-                        }
-                      double atmp = 0;
-                      for (octave_idx_type i = 0; i < j+1; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
                     }
-                  rcond = 1. / ainvnorm / anorm;
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-            }
-
-        triangular_error:
-          if (err != 0)
-            {
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
-            }
-
-          volatile double rcond_plus_one = rcond + 1.0;
-
-          if (rcond_plus_one == 1.0 || xisnan (rcond))
-            {
-              err = -2;
-
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
+              rcond = 1. / ainvnorm / anorm;
             }
         }
-      else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+    triangular_error:
+      if (err != 0)
+        {
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
+
+      volatile double rcond_plus_one = rcond + 1.0;
+
+      if (rcond_plus_one == 1.0 || xisnan (rcond))
+        {
+          err = -2;
+
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
     }
 
   return retval;
@@ -2722,7 +2684,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = Matrix (nc, b.cols (), 0.0);
   else
     {
@@ -2730,39 +2693,87 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower)
+      if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      double anorm = 0.;
+      double ainvnorm = 0.;
+      octave_idx_type b_nc = b.cols ();
+      rcond = 1.;
+
+      if (calc_cond)
+        {
+          // Calculate the 1-norm of matrix for rcond calculation
+          for (octave_idx_type j = 0; j < nc; j++)
+            {
+              double atmp = 0.;
+              for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                atmp += fabs (data (i));
+              if (atmp > anorm)
+                anorm = atmp;
+            }
+        }
+
+      if (typ == MatrixType::Permuted_Lower)
         {
-          double anorm = 0.;
-          double ainvnorm = 0.;
-          octave_idx_type b_nc = b.cols ();
-          rcond = 1.;
+          retval.resize (nc, b_nc);
+          OCTAVE_LOCAL_BUFFER (double, work, nm);
+          octave_idx_type *perm = mattype.triangular_perm ();
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              if (nc > nr)
+                for (octave_idx_type i = 0; i < nm; i++)
+                  work[i] = 0.;
+              for (octave_idx_type i = 0; i < nr; i++)
+                work[perm[i]] = b(i,j);
+
+              for (octave_idx_type k = 0; k < nc; k++)
+                {
+                  if (work[k] != 0.)
+                    {
+                      octave_idx_type minr = nr;
+                      octave_idx_type mini = 0;
+
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                        if (perm[ridx (i)] < minr)
+                          {
+                            minr = perm[ridx (i)];
+                            mini = i;
+                          }
+
+                      if (minr != k || data (mini) == 0)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      double tmp = work[k] / data (mini);
+                      work[k] = tmp;
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                        {
+                          if (i == mini)
+                            continue;
+
+                          octave_idx_type iidx = perm[ridx (i)];
+                          work[iidx] = work[iidx] - tmp * data (i);
+                        }
+                    }
+                }
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                retval(i, j) = work[i];
+            }
 
           if (calc_cond)
             {
-              // Calculate the 1-norm of matrix for rcond calculation
-              for (octave_idx_type j = 0; j < nc; j++)
+              // Calculation of 1-norm of inv(*this)
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  double atmp = 0.;
-                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
-                    atmp += fabs (data (i));
-                  if (atmp > anorm)
-                    anorm = atmp;
-                }
-            }
-
-          if (typ == MatrixType::Permuted_Lower)
-            {
-              retval.resize (nc, b_nc);
-              OCTAVE_LOCAL_BUFFER (double, work, nm);
-              octave_idx_type *perm = mattype.triangular_perm ();
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
-                {
-                  if (nc > nr)
-                    for (octave_idx_type i = 0; i < nm; i++)
-                      work[i] = 0.;
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    work[perm[i]] = b(i,j);
+                  work[j] = 1.;
 
                   for (octave_idx_type k = 0; k < nc; k++)
                     {
@@ -2771,22 +2782,18 @@
                           octave_idx_type minr = nr;
                           octave_idx_type mini = 0;
 
-                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                          for (octave_idx_type i = cidx (k);
+                               i < cidx (k+1); i++)
                             if (perm[ridx (i)] < minr)
                               {
                                 minr = perm[ridx (i)];
                                 mini = i;
                               }
 
-                          if (minr != k || data (mini) == 0)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
                           double tmp = work[k] / data (mini);
                           work[k] = tmp;
-                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                          for (octave_idx_type i = cidx (k);
+                               i < cidx (k+1); i++)
                             {
                               if (i == mini)
                                 continue;
@@ -2797,82 +2804,69 @@
                         }
                     }
 
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    retval(i, j) = work[i];
+                  double atmp = 0;
+                  for (octave_idx_type i = j; i < nc; i++)
+                    {
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
+                    }
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-
-              if (calc_cond)
+              rcond = 1. / ainvnorm / anorm;
+            }
+        }
+      else
+        {
+          OCTAVE_LOCAL_BUFFER (double, work, nm);
+          retval.resize (nc, b_nc, 0.);
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nr; i++)
+                work[i] = b(i,j);
+              for (octave_idx_type i = nr; i < nc; i++)
+                work[i] = 0.;
+              for (octave_idx_type k = 0; k < nc; k++)
                 {
-                  // Calculation of 1-norm of inv(*this)
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  if (work[k] != 0.)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = 0; k < nc; k++)
+                      if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
                         {
-                          if (work[k] != 0.)
-                            {
-                              octave_idx_type minr = nr;
-                              octave_idx_type mini = 0;
-
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
-                                if (perm[ridx (i)] < minr)
-                                  {
-                                    minr = perm[ridx (i)];
-                                    mini = i;
-                                  }
-
-                              double tmp = work[k] / data (mini);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
-                                {
-                                  if (i == mini)
-                                    continue;
-
-                                  octave_idx_type iidx = perm[ridx (i)];
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      double tmp = work[k] / data (cidx (k));
+                      work[k] = tmp;
+                      for (octave_idx_type i = cidx (k)+1;
+                           i < cidx (k+1); i++)
+                        {
+                          octave_idx_type iidx = ridx (i);
+                          work[iidx] = work[iidx] - tmp * data (i);
                         }
-
-                      double atmp = 0;
-                      for (octave_idx_type i = j; i < nc; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
                     }
-                  rcond = 1. / ainvnorm / anorm;
                 }
-            }
-          else
-            {
-              OCTAVE_LOCAL_BUFFER (double, work, nm);
-              retval.resize (nc, b_nc, 0.);
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                retval.xelem (i, j) = work[i];
+            }
+
+          if (calc_cond)
+            {
+              // Calculation of 1-norm of inv(*this)
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    work[i] = b(i,j);
-                  for (octave_idx_type i = nr; i < nc; i++)
-                    work[i] = 0.;
-                  for (octave_idx_type k = 0; k < nc; k++)
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k < nc; k++)
                     {
+
                       if (work[k] != 0.)
                         {
-                          if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
                           double tmp = work[k] / data (cidx (k));
                           work[k] = tmp;
                           for (octave_idx_type i = cidx (k)+1;
@@ -2883,78 +2877,45 @@
                             }
                         }
                     }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    retval.xelem (i, j) = work[i];
-                }
-
-              if (calc_cond)
-                {
-                  // Calculation of 1-norm of inv(*this)
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  double atmp = 0;
+                  for (octave_idx_type i = j; i < nc; i++)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k < nc; k++)
-                        {
-
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (k));
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
-                        }
-                      double atmp = 0;
-                      for (octave_idx_type i = j; i < nc; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
                     }
-                  rcond = 1. / ainvnorm / anorm;
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-            }
-
-        triangular_error:
-          if (err != 0)
-            {
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
-            }
-
-          volatile double rcond_plus_one = rcond + 1.0;
-
-          if (rcond_plus_one == 1.0 || xisnan (rcond))
-            {
-              err = -2;
-
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
+              rcond = 1. / ainvnorm / anorm;
             }
         }
-      else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+    triangular_error:
+      if (err != 0)
+        {
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
+
+      volatile double rcond_plus_one = rcond + 1.0;
+
+      if (rcond_plus_one == 1.0 || xisnan (rcond))
+        {
+          err = -2;
+
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
     }
 
   return retval;
@@ -2976,7 +2937,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = SparseMatrix (nc, b.cols ());
   else
     {
@@ -2984,43 +2946,113 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower)
+      if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      double anorm = 0.;
+      double ainvnorm = 0.;
+      rcond = 1.;
+
+      if (calc_cond)
+        {
+          // Calculate the 1-norm of matrix for rcond calculation
+          for (octave_idx_type j = 0; j < nc; j++)
+            {
+              double atmp = 0.;
+              for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                atmp += fabs (data (i));
+              if (atmp > anorm)
+                anorm = atmp;
+            }
+        }
+
+      octave_idx_type b_nc = b.cols ();
+      octave_idx_type b_nz = b.nnz ();
+      retval = SparseMatrix (nc, b_nc, b_nz);
+      retval.xcidx (0) = 0;
+      octave_idx_type ii = 0;
+      octave_idx_type x_nz = b_nz;
+
+      if (typ == MatrixType::Permuted_Lower)
         {
-          double anorm = 0.;
-          double ainvnorm = 0.;
-          rcond = 1.;
+          OCTAVE_LOCAL_BUFFER (double, work, nm);
+          octave_idx_type *perm = mattype.triangular_perm ();
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+              for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                work[perm[b.ridx (i)]] = b.data (i);
+
+              for (octave_idx_type k = 0; k < nc; k++)
+                {
+                  if (work[k] != 0.)
+                    {
+                      octave_idx_type minr = nr;
+                      octave_idx_type mini = 0;
+
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                        if (perm[ridx (i)] < minr)
+                          {
+                            minr = perm[ridx (i)];
+                            mini = i;
+                          }
+
+                      if (minr != k || data (mini) == 0)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      double tmp = work[k] / data (mini);
+                      work[k] = tmp;
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                        {
+                          if (i == mini)
+                            continue;
+
+                          octave_idx_type iidx = perm[ridx (i)];
+                          work[iidx] = work[iidx] - tmp * data (i);
+                        }
+                    }
+                }
+
+              // Count nonzeros in work vector and adjust space in
+              // retval if needed
+              octave_idx_type new_nnz = 0;
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (work[i] != 0.)
+                  new_nnz++;
+
+              if (ii + new_nnz > x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                  retval.change_capacity (sz);
+                  x_nz = sz;
+                }
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (work[i] != 0.)
+                  {
+                    retval.xridx (ii) = i;
+                    retval.xdata (ii++) = work[i];
+                  }
+              retval.xcidx (j+1) = ii;
+            }
+
+          retval.maybe_compress ();
 
           if (calc_cond)
             {
-              // Calculate the 1-norm of matrix for rcond calculation
-              for (octave_idx_type j = 0; j < nc; j++)
+              // Calculation of 1-norm of inv(*this)
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  double atmp = 0.;
-                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
-                    atmp += fabs (data (i));
-                  if (atmp > anorm)
-                    anorm = atmp;
-                }
-            }
-
-          octave_idx_type b_nc = b.cols ();
-          octave_idx_type b_nz = b.nnz ();
-          retval = SparseMatrix (nc, b_nc, b_nz);
-          retval.xcidx (0) = 0;
-          octave_idx_type ii = 0;
-          octave_idx_type x_nz = b_nz;
-
-          if (typ == MatrixType::Permuted_Lower)
-            {
-              OCTAVE_LOCAL_BUFFER (double, work, nm);
-              octave_idx_type *perm = mattype.triangular_perm ();
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
-                {
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
-                    work[perm[b.ridx (i)]] = b.data (i);
+                  work[j] = 1.;
 
                   for (octave_idx_type k = 0; k < nc; k++)
                     {
@@ -3029,22 +3061,18 @@
                           octave_idx_type minr = nr;
                           octave_idx_type mini = 0;
 
-                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                          for (octave_idx_type i = cidx (k);
+                               i < cidx (k+1); i++)
                             if (perm[ridx (i)] < minr)
                               {
                                 minr = perm[ridx (i)];
                                 mini = i;
                               }
 
-                          if (minr != k || data (mini) == 0)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
                           double tmp = work[k] / data (mini);
                           work[k] = tmp;
-                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                          for (octave_idx_type i = cidx (k);
+                               i < cidx (k+1); i++)
                             {
                               if (i == mini)
                                 continue;
@@ -3055,207 +3083,139 @@
                         }
                     }
 
-                  // Count nonzeros in work vector and adjust space in
-                  // retval if needed
-                  octave_idx_type new_nnz = 0;
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (work[i] != 0.)
-                      new_nnz++;
-
-                  if (ii + new_nnz > x_nz)
+                  double atmp = 0;
+                  for (octave_idx_type i = j; i < nr; i++)
                     {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
-                      retval.change_capacity (sz);
-                      x_nz = sz;
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
                     }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (work[i] != 0.)
-                      {
-                        retval.xridx (ii) = i;
-                        retval.xdata (ii++) = work[i];
-                      }
-                  retval.xcidx (j+1) = ii;
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-
-              retval.maybe_compress ();
-
-              if (calc_cond)
+              rcond = 1. / ainvnorm / anorm;
+            }
+        }
+      else
+        {
+          OCTAVE_LOCAL_BUFFER (double, work, nm);
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+              for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                work[b.ridx (i)] = b.data (i);
+
+              for (octave_idx_type k = 0; k < nc; k++)
                 {
-                  // Calculation of 1-norm of inv(*this)
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  if (work[k] != 0.)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = 0; k < nc; k++)
+                      if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      double tmp = work[k] / data (cidx (k));
+                      work[k] = tmp;
+                      for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
                         {
-                          if (work[k] != 0.)
-                            {
-                              octave_idx_type minr = nr;
-                              octave_idx_type mini = 0;
-
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
-                                if (perm[ridx (i)] < minr)
-                                  {
-                                    minr = perm[ridx (i)];
-                                    mini = i;
-                                  }
-
-                              double tmp = work[k] / data (mini);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
-                                {
-                                  if (i == mini)
-                                    continue;
-
-                                  octave_idx_type iidx = perm[ridx (i)];
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
+                          octave_idx_type iidx = ridx (i);
+                          work[iidx] = work[iidx] - tmp * data (i);
                         }
-
-                      double atmp = 0;
-                      for (octave_idx_type i = j; i < nr; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
                     }
-                  rcond = 1. / ainvnorm / anorm;
+                }
+
+              // Count nonzeros in work vector and adjust space in
+              // retval if needed
+              octave_idx_type new_nnz = 0;
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (work[i] != 0.)
+                  new_nnz++;
+
+              if (ii + new_nnz > x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                  retval.change_capacity (sz);
+                  x_nz = sz;
                 }
-            }
-          else
-            {
-              OCTAVE_LOCAL_BUFFER (double, work, nm);
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (work[i] != 0.)
+                  {
+                    retval.xridx (ii) = i;
+                    retval.xdata (ii++) = work[i];
+                  }
+              retval.xcidx (j+1) = ii;
+            }
+
+          retval.maybe_compress ();
+
+          if (calc_cond)
+            {
+              // Calculation of 1-norm of inv(*this)
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
-                    work[b.ridx (i)] = b.data (i);
-
-                  for (octave_idx_type k = 0; k < nc; k++)
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k < nc; k++)
                     {
+
                       if (work[k] != 0.)
                         {
-                          if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
                           double tmp = work[k] / data (cidx (k));
                           work[k] = tmp;
-                          for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                          for (octave_idx_type i = cidx (k)+1;
+                               i < cidx (k+1); i++)
                             {
                               octave_idx_type iidx = ridx (i);
                               work[iidx] = work[iidx] - tmp * data (i);
                             }
                         }
                     }
-
-                  // Count nonzeros in work vector and adjust space in
-                  // retval if needed
-                  octave_idx_type new_nnz = 0;
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (work[i] != 0.)
-                      new_nnz++;
-
-                  if (ii + new_nnz > x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
-                      retval.change_capacity (sz);
-                      x_nz = sz;
-                    }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (work[i] != 0.)
-                      {
-                        retval.xridx (ii) = i;
-                        retval.xdata (ii++) = work[i];
-                      }
-                  retval.xcidx (j+1) = ii;
-                }
-
-              retval.maybe_compress ();
-
-              if (calc_cond)
-                {
-                  // Calculation of 1-norm of inv(*this)
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  double atmp = 0;
+                  for (octave_idx_type i = j; i < nc; i++)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k < nc; k++)
-                        {
-
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (k));
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
-                        }
-                      double atmp = 0;
-                      for (octave_idx_type i = j; i < nc; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
                     }
-                  rcond = 1. / ainvnorm / anorm;
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-            }
-
-        triangular_error:
-          if (err != 0)
-            {
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
-            }
-
-          volatile double rcond_plus_one = rcond + 1.0;
-
-          if (rcond_plus_one == 1.0 || xisnan (rcond))
-            {
-              err = -2;
-
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
+              rcond = 1. / ainvnorm / anorm;
             }
         }
-      else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+    triangular_error:
+      if (err != 0)
+        {
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
+
+      volatile double rcond_plus_one = rcond + 1.0;
+
+      if (rcond_plus_one == 1.0 || xisnan (rcond))
+        {
+          err = -2;
+
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
     }
 
   return retval;
@@ -3277,7 +3237,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
   else
     {
@@ -3285,232 +3246,230 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower)
+      if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      double anorm = 0.;
+      double ainvnorm = 0.;
+      octave_idx_type b_nc = b.cols ();
+      rcond = 1.;
+
+      if (calc_cond)
+        {
+          // Calculate the 1-norm of matrix for rcond calculation
+          for (octave_idx_type j = 0; j < nc; j++)
+            {
+              double atmp = 0.;
+              for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                atmp += fabs (data (i));
+              if (atmp > anorm)
+                anorm = atmp;
+            }
+        }
+
+      if (typ == MatrixType::Permuted_Lower)
         {
-          double anorm = 0.;
-          double ainvnorm = 0.;
-          octave_idx_type b_nc = b.cols ();
-          rcond = 1.;
+          retval.resize (nc, b_nc);
+          OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+          octave_idx_type *perm = mattype.triangular_perm ();
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nm; i++)
+                cwork[i] = 0.;
+              for (octave_idx_type i = 0; i < nr; i++)
+                cwork[perm[i]] = b(i,j);
+
+              for (octave_idx_type k = 0; k < nc; k++)
+                {
+                  if (cwork[k] != 0.)
+                    {
+                      octave_idx_type minr = nr;
+                      octave_idx_type mini = 0;
+
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                        if (perm[ridx (i)] < minr)
+                          {
+                            minr = perm[ridx (i)];
+                            mini = i;
+                          }
+
+                      if (minr != k || data (mini) == 0)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      Complex tmp = cwork[k] / data (mini);
+                      cwork[k] = tmp;
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                        {
+                          if (i == mini)
+                            continue;
+
+                          octave_idx_type iidx = perm[ridx (i)];
+                          cwork[iidx] = cwork[iidx] - tmp * data (i);
+                        }
+                    }
+                }
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                retval(i, j) = cwork[i];
+            }
 
           if (calc_cond)
             {
-              // Calculate the 1-norm of matrix for rcond calculation
-              for (octave_idx_type j = 0; j < nc; j++)
+              // Calculation of 1-norm of inv(*this)
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  double atmp = 0.;
-                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
-                    atmp += fabs (data (i));
-                  if (atmp > anorm)
-                    anorm = atmp;
-                }
-            }
-
-          if (typ == MatrixType::Permuted_Lower)
-            {
-              retval.resize (nc, b_nc);
-              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
-              octave_idx_type *perm = mattype.triangular_perm ();
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
-                {
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    cwork[i] = 0.;
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    cwork[perm[i]] = b(i,j);
+                  work[j] = 1.;
 
                   for (octave_idx_type k = 0; k < nc; k++)
                     {
-                      if (cwork[k] != 0.)
+                      if (work[k] != 0.)
                         {
                           octave_idx_type minr = nr;
                           octave_idx_type mini = 0;
 
-                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                          for (octave_idx_type i = cidx (k);
+                               i < cidx (k+1); i++)
                             if (perm[ridx (i)] < minr)
                               {
                                 minr = perm[ridx (i)];
                                 mini = i;
                               }
 
-                          if (minr != k || data (mini) == 0)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
-                          Complex tmp = cwork[k] / data (mini);
-                          cwork[k] = tmp;
-                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                          double tmp = work[k] / data (mini);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k);
+                               i < cidx (k+1); i++)
                             {
                               if (i == mini)
                                 continue;
 
                               octave_idx_type iidx = perm[ridx (i)];
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
                             }
                         }
                     }
 
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    retval(i, j) = cwork[i];
+                  double atmp = 0;
+                  for (octave_idx_type i = j; i < nc; i++)
+                    {
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
+                    }
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-
-              if (calc_cond)
+              rcond = 1. / ainvnorm / anorm;
+            }
+        }
+      else
+        {
+          OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+          retval.resize (nc, b_nc, 0.);
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nr; i++)
+                cwork[i] = b(i,j);
+              for (octave_idx_type i = nr; i < nc; i++)
+                cwork[i] = 0.;
+
+              for (octave_idx_type k = 0; k < nc; k++)
                 {
-                  // Calculation of 1-norm of inv(*this)
-                  OCTAVE_LOCAL_BUFFER (double, work, nm);
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  if (cwork[k] != 0.)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = 0; k < nc; k++)
+                      if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
                         {
-                          if (work[k] != 0.)
-                            {
-                              octave_idx_type minr = nr;
-                              octave_idx_type mini = 0;
-
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
-                                if (perm[ridx (i)] < minr)
-                                  {
-                                    minr = perm[ridx (i)];
-                                    mini = i;
-                                  }
-
-                              double tmp = work[k] / data (mini);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
-                                {
-                                  if (i == mini)
-                                    continue;
-
-                                  octave_idx_type iidx = perm[ridx (i)];
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
+                          err = -2;
+                          goto triangular_error;
                         }
 
-                      double atmp = 0;
-                      for (octave_idx_type i = j; i < nc; i++)
+                      Complex tmp = cwork[k] / data (cidx (k));
+                      cwork[k] = tmp;
+                      for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
                         {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
+                          octave_idx_type iidx = ridx (i);
+                          cwork[iidx] = cwork[iidx] - tmp * data (i);
                         }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
                     }
-                  rcond = 1. / ainvnorm / anorm;
                 }
-            }
-          else
-            {
-              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
-              retval.resize (nc, b_nc, 0.);
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                retval.xelem (i, j) = cwork[i];
+            }
+
+          if (calc_cond)
+            {
+              // Calculation of 1-norm of inv(*this)
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  for (octave_idx_type i = 0; i < nr; i++)
-                    cwork[i] = b(i,j);
-                  for (octave_idx_type i = nr; i < nc; i++)
-                    cwork[i] = 0.;
-
-                  for (octave_idx_type k = 0; k < nc; k++)
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k < nc; k++)
                     {
-                      if (cwork[k] != 0.)
+
+                      if (work[k] != 0.)
                         {
-                          if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
-                          Complex tmp = cwork[k] / data (cidx (k));
-                          cwork[k] = tmp;
-                          for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                          double tmp = work[k] / data (cidx (k));
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k)+1;
+                               i < cidx (k+1); i++)
                             {
                               octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
                             }
                         }
                     }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    retval.xelem (i, j) = cwork[i];
-                }
-
-              if (calc_cond)
-                {
-                  // Calculation of 1-norm of inv(*this)
-                  OCTAVE_LOCAL_BUFFER (double, work, nm);
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  double atmp = 0;
+                  for (octave_idx_type i = j; i < nc; i++)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k < nc; k++)
-                        {
-
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (k));
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
-                        }
-                      double atmp = 0;
-                      for (octave_idx_type i = j; i < nc; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
                     }
-                  rcond = 1. / ainvnorm / anorm;
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-            }
-
-        triangular_error:
-          if (err != 0)
-            {
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
-            }
-
-          volatile double rcond_plus_one = rcond + 1.0;
-
-          if (rcond_plus_one == 1.0 || xisnan (rcond))
-            {
-              err = -2;
-
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
+              rcond = 1. / ainvnorm / anorm;
             }
         }
-      else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+    triangular_error:
+      if (err != 0)
+        {
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
+
+      volatile double rcond_plus_one = rcond + 1.0;
+
+      if (rcond_plus_one == 1.0 || xisnan (rcond))
+        {
+          err = -2;
+
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
     }
 
   return retval;
@@ -3532,7 +3491,8 @@
   if (nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || nc == 0 || b.cols () == 0)
+
+  if (nr == 0 || nc == 0 || b.cols () == 0)
     retval = SparseComplexMatrix (nc, b.cols ());
   else
     {
@@ -3540,280 +3500,278 @@
       int typ = mattype.type ();
       mattype.info ();
 
-      if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower)
+      if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
+        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+      double anorm = 0.;
+      double ainvnorm = 0.;
+      rcond = 1.;
+
+      if (calc_cond)
+        {
+          // Calculate the 1-norm of matrix for rcond calculation
+          for (octave_idx_type j = 0; j < nc; j++)
+            {
+              double atmp = 0.;
+              for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
+                atmp += fabs (data (i));
+              if (atmp > anorm)
+                anorm = atmp;
+            }
+        }
+
+      octave_idx_type b_nc = b.cols ();
+      octave_idx_type b_nz = b.nnz ();
+      retval = SparseComplexMatrix (nc, b_nc, b_nz);
+      retval.xcidx (0) = 0;
+      octave_idx_type ii = 0;
+      octave_idx_type x_nz = b_nz;
+
+      if (typ == MatrixType::Permuted_Lower)
         {
-          double anorm = 0.;
-          double ainvnorm = 0.;
-          rcond = 1.;
+          OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+          octave_idx_type *perm = mattype.triangular_perm ();
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nm; i++)
+                cwork[i] = 0.;
+              for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                cwork[perm[b.ridx (i)]] = b.data (i);
+
+              for (octave_idx_type k = 0; k < nc; k++)
+                {
+                  if (cwork[k] != 0.)
+                    {
+                      octave_idx_type minr = nr;
+                      octave_idx_type mini = 0;
+
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                        if (perm[ridx (i)] < minr)
+                          {
+                            minr = perm[ridx (i)];
+                            mini = i;
+                          }
+
+                      if (minr != k || data (mini) == 0)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      Complex tmp = cwork[k] / data (mini);
+                      cwork[k] = tmp;
+                      for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                        {
+                          if (i == mini)
+                            continue;
+
+                          octave_idx_type iidx = perm[ridx (i)];
+                          cwork[iidx] = cwork[iidx] - tmp * data (i);
+                        }
+                    }
+                }
+
+              // Count nonzeros in work vector and adjust space in
+              // retval if needed
+              octave_idx_type new_nnz = 0;
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (cwork[i] != 0.)
+                  new_nnz++;
+
+              if (ii + new_nnz > x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                  retval.change_capacity (sz);
+                  x_nz = sz;
+                }
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (cwork[i] != 0.)
+                  {
+                    retval.xridx (ii) = i;
+                    retval.xdata (ii++) = cwork[i];
+                  }
+              retval.xcidx (j+1) = ii;
+            }
+
+          retval.maybe_compress ();
 
           if (calc_cond)
             {
-              // Calculate the 1-norm of matrix for rcond calculation
-              for (octave_idx_type j = 0; j < nc; j++)
+              // Calculation of 1-norm of inv(*this)
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
                 {
-                  double atmp = 0.;
-                  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
-                    atmp += fabs (data (i));
-                  if (atmp > anorm)
-                    anorm = atmp;
-                }
-            }
-
-          octave_idx_type b_nc = b.cols ();
-          octave_idx_type b_nz = b.nnz ();
-          retval = SparseComplexMatrix (nc, b_nc, b_nz);
-          retval.xcidx (0) = 0;
-          octave_idx_type ii = 0;
-          octave_idx_type x_nz = b_nz;
-
-          if (typ == MatrixType::Permuted_Lower)
-            {
-              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
-              octave_idx_type *perm = mattype.triangular_perm ();
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
-                {
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    cwork[i] = 0.;
-                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
-                    cwork[perm[b.ridx (i)]] = b.data (i);
+                  work[j] = 1.;
 
                   for (octave_idx_type k = 0; k < nc; k++)
                     {
-                      if (cwork[k] != 0.)
+                      if (work[k] != 0.)
                         {
                           octave_idx_type minr = nr;
                           octave_idx_type mini = 0;
 
-                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                          for (octave_idx_type i = cidx (k);
+                               i < cidx (k+1); i++)
                             if (perm[ridx (i)] < minr)
                               {
                                 minr = perm[ridx (i)];
                                 mini = i;
                               }
 
-                          if (minr != k || data (mini) == 0)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
-                          Complex tmp = cwork[k] / data (mini);
-                          cwork[k] = tmp;
-                          for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
+                          double tmp = work[k] / data (mini);
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k);
+                               i < cidx (k+1); i++)
                             {
                               if (i == mini)
                                 continue;
 
                               octave_idx_type iidx = perm[ridx (i)];
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
                             }
                         }
                     }
 
-                  // Count nonzeros in work vector and adjust space in
-                  // retval if needed
-                  octave_idx_type new_nnz = 0;
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (cwork[i] != 0.)
-                      new_nnz++;
-
-                  if (ii + new_nnz > x_nz)
+                  double atmp = 0;
+                  for (octave_idx_type i = j; i < nc; i++)
                     {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
-                      retval.change_capacity (sz);
-                      x_nz = sz;
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
                     }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (cwork[i] != 0.)
-                      {
-                        retval.xridx (ii) = i;
-                        retval.xdata (ii++) = cwork[i];
-                      }
-                  retval.xcidx (j+1) = ii;
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-
-              retval.maybe_compress ();
-
-              if (calc_cond)
+              rcond = 1. / ainvnorm / anorm;
+            }
+        }
+      else
+        {
+          OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
+
+          for (octave_idx_type j = 0; j < b_nc; j++)
+            {
+              for (octave_idx_type i = 0; i < nm; i++)
+                cwork[i] = 0.;
+              for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
+                cwork[b.ridx (i)] = b.data (i);
+
+              for (octave_idx_type k = 0; k < nc; k++)
                 {
-                  // Calculation of 1-norm of inv(*this)
-                  OCTAVE_LOCAL_BUFFER (double, work, nm);
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  if (cwork[k] != 0.)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = 0; k < nc; k++)
+                      if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
+                        {
+                          err = -2;
+                          goto triangular_error;
+                        }
+
+                      Complex tmp = cwork[k] / data (cidx (k));
+                      cwork[k] = tmp;
+                      for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
                         {
-                          if (work[k] != 0.)
-                            {
-                              octave_idx_type minr = nr;
-                              octave_idx_type mini = 0;
-
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
-                                if (perm[ridx (i)] < minr)
-                                  {
-                                    minr = perm[ridx (i)];
-                                    mini = i;
-                                  }
-
-                              double tmp = work[k] / data (mini);
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k);
-                                   i < cidx (k+1); i++)
-                                {
-                                  if (i == mini)
-                                    continue;
-
-                                  octave_idx_type iidx = perm[ridx (i)];
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
+                          octave_idx_type iidx = ridx (i);
+                          cwork[iidx] = cwork[iidx] - tmp * data (i);
                         }
-
-                      double atmp = 0;
-                      for (octave_idx_type i = j; i < nc; i++)
+                    }
+                }
+
+              // Count nonzeros in work vector and adjust space in
+              // retval if needed
+              octave_idx_type new_nnz = 0;
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (cwork[i] != 0.)
+                  new_nnz++;
+
+              if (ii + new_nnz > x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
+                  retval.change_capacity (sz);
+                  x_nz = sz;
+                }
+
+              for (octave_idx_type i = 0; i < nc; i++)
+                if (cwork[i] != 0.)
+                  {
+                    retval.xridx (ii) = i;
+                    retval.xdata (ii++) = cwork[i];
+                  }
+              retval.xcidx (j+1) = ii;
+            }
+
+          retval.maybe_compress ();
+
+          if (calc_cond)
+            {
+              // Calculation of 1-norm of inv(*this)
+              OCTAVE_LOCAL_BUFFER (double, work, nm);
+              for (octave_idx_type i = 0; i < nm; i++)
+                work[i] = 0.;
+
+              for (octave_idx_type j = 0; j < nr; j++)
+                {
+                  work[j] = 1.;
+
+                  for (octave_idx_type k = j; k < nc; k++)
+                    {
+
+                      if (work[k] != 0.)
                         {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
-                    }
-                  rcond = 1. / ainvnorm / anorm;
-                }
-            }
-          else
-            {
-              OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
-
-              for (octave_idx_type j = 0; j < b_nc; j++)
-                {
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    cwork[i] = 0.;
-                  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
-                    cwork[b.ridx (i)] = b.data (i);
-
-                  for (octave_idx_type k = 0; k < nc; k++)
-                    {
-                      if (cwork[k] != 0.)
-                        {
-                          if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
-                            {
-                              err = -2;
-                              goto triangular_error;
-                            }
-
-                          Complex tmp = cwork[k] / data (cidx (k));
-                          cwork[k] = tmp;
-                          for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
+                          double tmp = work[k] / data (cidx (k));
+                          work[k] = tmp;
+                          for (octave_idx_type i = cidx (k)+1;
+                               i < cidx (k+1); i++)
                             {
                               octave_idx_type iidx = ridx (i);
-                              cwork[iidx] = cwork[iidx] - tmp * data (i);
+                              work[iidx] = work[iidx] - tmp * data (i);
                             }
                         }
                     }
-
-                  // Count nonzeros in work vector and adjust space in
-                  // retval if needed
-                  octave_idx_type new_nnz = 0;
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (cwork[i] != 0.)
-                      new_nnz++;
-
-                  if (ii + new_nnz > x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
-                      retval.change_capacity (sz);
-                      x_nz = sz;
-                    }
-
-                  for (octave_idx_type i = 0; i < nc; i++)
-                    if (cwork[i] != 0.)
-                      {
-                        retval.xridx (ii) = i;
-                        retval.xdata (ii++) = cwork[i];
-                      }
-                  retval.xcidx (j+1) = ii;
-                }
-
-              retval.maybe_compress ();
-
-              if (calc_cond)
-                {
-                  // Calculation of 1-norm of inv(*this)
-                  OCTAVE_LOCAL_BUFFER (double, work, nm);
-                  for (octave_idx_type i = 0; i < nm; i++)
-                    work[i] = 0.;
-
-                  for (octave_idx_type j = 0; j < nr; j++)
+                  double atmp = 0;
+                  for (octave_idx_type i = j; i < nc; i++)
                     {
-                      work[j] = 1.;
-
-                      for (octave_idx_type k = j; k < nc; k++)
-                        {
-
-                          if (work[k] != 0.)
-                            {
-                              double tmp = work[k] / data (cidx (k));
-                              work[k] = tmp;
-                              for (octave_idx_type i = cidx (k)+1;
-                                   i < cidx (k+1); i++)
-                                {
-                                  octave_idx_type iidx = ridx (i);
-                                  work[iidx] = work[iidx] - tmp * data (i);
-                                }
-                            }
-                        }
-                      double atmp = 0;
-                      for (octave_idx_type i = j; i < nc; i++)
-                        {
-                          atmp += fabs (work[i]);
-                          work[i] = 0.;
-                        }
-                      if (atmp > ainvnorm)
-                        ainvnorm = atmp;
+                      atmp += fabs (work[i]);
+                      work[i] = 0.;
                     }
-                  rcond = 1. / ainvnorm / anorm;
+                  if (atmp > ainvnorm)
+                    ainvnorm = atmp;
                 }
-            }
-
-        triangular_error:
-          if (err != 0)
-            {
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
-            }
-
-          volatile double rcond_plus_one = rcond + 1.0;
-
-          if (rcond_plus_one == 1.0 || xisnan (rcond))
-            {
-              err = -2;
-
-              if (sing_handler)
-                {
-                  sing_handler (rcond);
-                  mattype.mark_as_rectangular ();
-                }
-              else
-                warn_singular_matrix (rcond);
+              rcond = 1. / ainvnorm / anorm;
             }
         }
-      else
-        (*current_liboctave_error_handler) ("incorrect matrix type");
+
+    triangular_error:
+      if (err != 0)
+        {
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
+
+      volatile double rcond_plus_one = rcond + 1.0;
+
+      if (rcond_plus_one == 1.0 || xisnan (rcond))
+        {
+          err = -2;
+
+          if (sing_handler)
+            {
+              sing_handler (rcond);
+              mattype.mark_as_rectangular ();
+            }
+          else
+            warn_singular_matrix (rcond);
+        }
     }
 
   return retval;
@@ -3834,7 +3792,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = Matrix (nc, b.cols (), 0.0);
   else if (calc_cond)
     (*current_liboctave_error_handler)
@@ -3983,7 +3942,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = SparseMatrix (nc, b.cols ());
   else if (calc_cond)
     (*current_liboctave_error_handler)
@@ -4128,7 +4088,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
   else if (calc_cond)
     (*current_liboctave_error_handler)
@@ -4279,7 +4240,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = SparseComplexMatrix (nc, b.cols ());
   else if (calc_cond)
     (*current_liboctave_error_handler)
@@ -4384,6 +4346,7 @@
 
                   if (err != 0)
                     {
+                      // FIXME: Should this be a warning?
                       (*current_liboctave_error_handler)
                         ("SparseMatrix::solve solve failed");
 
@@ -4399,6 +4362,7 @@
 
                   if (err != 0)
                     {
+                      // FIXME: Should this be a warning?
                       (*current_liboctave_error_handler)
                         ("SparseMatrix::solve solve failed");
 
@@ -4457,7 +4421,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = Matrix (nc, b.cols (), 0.0);
   else
     {
@@ -4559,6 +4524,7 @@
 
                   if (err != 0)
                     {
+                      // FIXME: Should this be a warning?
                       (*current_liboctave_error_handler)
                         ("SparseMatrix::solve solve failed");
                       err = -1;
@@ -4700,7 +4666,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = SparseMatrix (nc, b.cols ());
   else
     {
@@ -4812,6 +4779,7 @@
 
                       if (err != 0)
                         {
+                          // FIXME: Should this be a warning?
                           (*current_liboctave_error_handler)
                             ("SparseMatrix::solve solve failed");
                           err = -1;
@@ -5012,7 +4980,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
   else
     {
@@ -5127,6 +5096,7 @@
 
                       if (err != 0)
                         {
+                          // FIXME: Should this be a warning?
                           (*current_liboctave_error_handler)
                             ("SparseMatrix::solve solve failed");
                           err = -1;
@@ -5141,6 +5111,7 @@
 
                       if (err != 0)
                         {
+                          // FIXME: Should this be a warning?
                           (*current_liboctave_error_handler)
                             ("SparseMatrix::solve solve failed");
                           err = -1;
@@ -5305,7 +5276,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = SparseComplexMatrix (nc, b.cols ());
   else
     {
@@ -5426,6 +5398,7 @@
 
                       if (err != 0)
                         {
+                          // FIXME: Should this be a warning?
                           (*current_liboctave_error_handler)
                             ("SparseMatrix::solve solve failed");
                           err = -1;
@@ -5440,6 +5413,7 @@
 
                       if (err != 0)
                         {
+                          // FIXME: Should this be a warning?
                           (*current_liboctave_error_handler)
                             ("SparseMatrix::solve solve failed");
 
@@ -5694,14 +5668,15 @@
 
   if (status < 0)
     {
-      (*current_liboctave_error_handler)
-        ("SparseMatrix::solve symbolic factorization failed");
-      err = -1;
-
       UMFPACK_DNAME (report_status) (control, status);
       UMFPACK_DNAME (report_info) (control, info);
 
       UMFPACK_DNAME (free_symbolic) (&Symbolic);
+
+      // FIXME: Should this be a warning?
+      (*current_liboctave_error_handler)
+        ("SparseMatrix::solve symbolic factorization failed");
+      err = -1;
     }
   else
     {
@@ -5731,12 +5706,13 @@
         }
       else if (status < 0)
         {
+          UMFPACK_DNAME (report_status) (control, status);
+          UMFPACK_DNAME (report_info) (control, info);
+
+          // FIXME: Should this be a warning?
           (*current_liboctave_error_handler)
             ("SparseMatrix::solve numeric factorization failed");
 
-          UMFPACK_DNAME (report_status) (control, status);
-          UMFPACK_DNAME (report_info) (control, info);
-
           err = -1;
         }
       else
@@ -5772,7 +5748,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = Matrix (nc, b.cols (), 0.0);
   else
     {
@@ -5945,13 +5922,13 @@
                                                   info);
                   if (status < 0)
                     {
+                      UMFPACK_DNAME (report_status) (control, status);
+
+                      // FIXME: Should this be a warning?
                       (*current_liboctave_error_handler)
                         ("SparseMatrix::solve solve failed");
 
-                      UMFPACK_DNAME (report_status) (control, status);
-
                       err = -1;
-
                       break;
                     }
                 }
@@ -5991,7 +5968,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = SparseMatrix (nc, b.cols ());
   else
     {
@@ -6188,13 +6166,13 @@
                                                   control, info);
                   if (status < 0)
                     {
+                      UMFPACK_DNAME (report_status) (control, status);
+
+                      // FIXME: Should this be a warning?
                       (*current_liboctave_error_handler)
                         ("SparseMatrix::solve solve failed");
 
-                      UMFPACK_DNAME (report_status) (control, status);
-
                       err = -1;
-
                       break;
                     }
 
@@ -6255,7 +6233,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
   else
     {
@@ -6442,13 +6421,13 @@
 
                   if (status < 0 || status2 < 0)
                     {
+                      UMFPACK_DNAME (report_status) (control, status);
+
+                      // FIXME: Should this be a warning?
                       (*current_liboctave_error_handler)
                         ("SparseMatrix::solve solve failed");
 
-                      UMFPACK_DNAME (report_status) (control, status);
-
                       err = -1;
-
                       break;
                     }
 
@@ -6491,7 +6470,8 @@
   if (nr != nc || nr != b.rows ())
     (*current_liboctave_error_handler)
       ("matrix dimension mismatch solution of linear equations");
-  else if (nr == 0 || b.cols () == 0)
+
+  if (nr == 0 || b.cols () == 0)
     retval = SparseComplexMatrix (nc, b.cols ());
   else
     {
@@ -6699,13 +6679,13 @@
 
                   if (status < 0 || status2 < 0)
                     {
+                      UMFPACK_DNAME (report_status) (control, status);
+
+                      // FIXME: Should this be a warning?
                       (*current_liboctave_error_handler)
                         ("SparseMatrix::solve solve failed");
 
-                      UMFPACK_DNAME (report_status) (control, status);
-
                       err = -1;
-
                       break;
                     }
 
@@ -6799,10 +6779,7 @@
   else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
     retval = fsolve (mattype, b, err, rcond, sing_handler, true);
   else if (typ != MatrixType::Rectangular)
-    {
-      (*current_liboctave_error_handler) ("unknown matrix type");
-      return Matrix ();
-    }
+    (*current_liboctave_error_handler) ("unknown matrix type");
 
   // Rectangular or one of the above solvers flags a singular matrix
   if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
@@ -6867,10 +6844,7 @@
   else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
     retval = fsolve (mattype, b, err, rcond, sing_handler, true);
   else if (typ != MatrixType::Rectangular)
-    {
-      (*current_liboctave_error_handler) ("unknown matrix type");
-      return SparseMatrix ();
-    }
+    (*current_liboctave_error_handler) ("unknown matrix type");
 
   if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
     {
@@ -6935,10 +6909,7 @@
   else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
     retval = fsolve (mattype, b, err, rcond, sing_handler, true);
   else if (typ != MatrixType::Rectangular)
-    {
-      (*current_liboctave_error_handler) ("unknown matrix type");
-      return ComplexMatrix ();
-    }
+    (*current_liboctave_error_handler) ("unknown matrix type");
 
   if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
     {
@@ -7003,10 +6974,7 @@
   else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
     retval = fsolve (mattype, b, err, rcond, sing_handler, true);
   else if (typ != MatrixType::Rectangular)
-    {
-      (*current_liboctave_error_handler) ("unknown matrix type");
-      return SparseComplexMatrix ();
-    }
+    (*current_liboctave_error_handler) ("unknown matrix type");
 
   if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
     {