diff liboctave/dSparse.cc @ 5876:565d0cd4d9d0

[project @ 2006-07-01 19:42:06 by dbateman]
author dbateman
date Sat, 01 Jul 2006 19:42:07 +0000
parents 6b9cec830d72
children 50d43cdbec80
line wrap: on
line diff
--- a/liboctave/dSparse.cc
+++ b/liboctave/dSparse.cc
@@ -827,6 +827,14 @@
 		      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;
@@ -847,11 +855,17 @@
 				 (colXp<cx) && (colUp<nz));
 
 		      // get A(m,m)
-		      colUp = cidx(j+1) - 1;
+		      if (typ == MatrixType::Upper)
+			colUp = cidx(j+1) - 1;
+		      else
+			colUp = cidx(j) - 1;
 		      double pivot = data(colUp);
-		      if (pivot == 0.) 
-			(*current_liboctave_error_handler) 
-			  ("division by zero");
+		      if (pivot == 0. || colUp != j) 
+			{
+			  (*current_liboctave_error_handler) 
+			    ("division by zero");
+			  goto inverse_singular;
+			}
 
 		      if (v != 0.)
 			{
@@ -868,10 +882,17 @@
 		    }
 
 		  // get A(m,m)
-		  octave_idx_type colUp = cidx(i+1) - 1;
+		  octave_idx_type colUp;
+		  if (typ == MatrixType::Upper)
+		    colUp = cidx(i+1) - 1;
+		  else
+		    colUp = cidx(i) - 1;
 		  double pivot = data(colUp);
-		  if (pivot == 0.) 
-		    (*current_liboctave_error_handler) ("division by zero");
+		  if (pivot == 0. || 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++)
@@ -929,20 +950,35 @@
 			}
 
 		      // get A(m,m)
-		      double pivot = data(cidx(jidx+1) - 1);
+		      double pivot;
+		      if (typ == MatrixType::Permuted_Upper)
+			pivot = data(cidx(jidx+1) - 1);
+		      else
+			pivot = data(cidx(jidx) - 1);
 		      if (pivot == 0.) 
-			(*current_liboctave_error_handler) 
-			  ("division by zero");
+			{
+			  (*current_liboctave_error_handler) 
+			    ("division by zero");
+			  goto inverse_singular;
+			}
 
 		      work[j] = v / pivot;
 		    }
 
 		  // get A(m,m)
-		  octave_idx_type colUp = cidx(perm[iidx]+1) - 1;
+		  octave_idx_type colUp;
+		  if (typ == MatrixType::Permuted_Upper)
+		    colUp = cidx(perm[iidx]+1) - 1;
+		  else
+		    colUp = cidx(perm[iidx]) - 1;		  
+
 		  double pivot = data(colUp);
-		  if (pivot == 0.) 
-		    (*current_liboctave_error_handler) 
-		      ("division by zero");
+		  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++)
@@ -993,6 +1029,9 @@
     }
 
   return retval;
+
+ inverse_singular:
+  return SparseMatrix();
 }
 
 SparseMatrix