Mercurial > hg > octave-thorsten

changeset 5876:565d0cd4d9d0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-07-01 19:42:06 by dbateman]
author dbateman
date Sat, 01 Jul 2006 19:42:07 +0000
parents f6ddc0ee2315
children 50d43cdbec80
files liboctave/CSparse.cc liboctave/ChangeLog liboctave/Sparse-op-defs.h liboctave/dSparse.cc
diffstat 4 files changed, 118 insertions(+), 30 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CSparse.cc
+++ b/liboctave/CSparse.cc
@@ -750,6 +750,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;
@@ -769,12 +777,19 @@
 			} while ((rpX<j) && (rpU<j) && 
 				 (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;
 		      Complex 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.)
 			{
@@ -791,10 +806,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;
 		  Complex 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++)
@@ -852,20 +874,35 @@
 			}
 
 		      // get A(m,m)
-		      Complex pivot = data(cidx(jidx+1) - 1);
+		      Complex 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;
-		  Complex pivot = data(colUp);
-		  if (pivot == 0.) 
-		    (*current_liboctave_error_handler) 
-		      ("division by zero");
+		  octave_idx_type colUp;
+		  if (typ == MatrixType::Permuted_Upper)
+		    colUp = cidx(perm[iidx]+1) - 1;
+		  else
+		    colUp = cidx(perm[iidx]) - 1;		  
+
+  		  Complex 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++)
@@ -916,6 +953,9 @@
     }
 
   return retval;
+
+ inverse_singular:
+  return SparseComplexMatrix();
 }
 
 SparseComplexMatrix
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,12 @@
+2006-07-01  David Bateman  <dbateman@free.fr>
+
+	* dSparse.cc (tinverse): Check for rows with no elements and zero
+	elements on the diagonal. Allow both Upper and Lower triangular
+	matrices to be treated.
+	* CSparse.cc (tinverse): ditto.
+	* Sparse-op-defs.h (SPARSE_SPARSE_MUL): Take into account 64-bit
+	constant assignment.
+	
 2006-06-30  John W. Eaton  <jwe@octave.org>
 
 	* lo-sysdep.cc (octave_chdir): Perform tilde expansion here.
--- a/liboctave/Sparse-op-defs.h
+++ b/liboctave/Sparse-op-defs.h
@@ -1548,7 +1548,7 @@
   else \
     { \
       OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr); \
-      RET_TYPE retval (nr, a_nc, 0); \
+      RET_TYPE retval (nr, a_nc, static_cast<octave_idx_type> (0)); \
       for (octave_idx_type i = 0; i < nr; i++) \
 	w[i] = 0; \
       retval.xcidx(0) = 0; \
--- 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