Mercurial > hg > octave-nkf

diff liboctave/CMatrix.cc @ 7076:0bade2dc44a1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-29 18:09:57 by jwe]
author jwe
date Mon, 29 Oct 2007 18:09:57 +0000
parents b48d486f641d
children 6d3e53a2f963
line wrap: on
line diff
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -2201,7 +2201,7 @@
   if (singular_fallback && mattype.type () == MatrixType::Rectangular)
     {
       octave_idx_type rank;
-      retval = lssolve (b, info, rank);
+      retval = lssolve (b, info, rank, rcond);
     }
 
   return retval;
@@ -2395,20 +2395,31 @@
 {
   octave_idx_type info;
   octave_idx_type rank;
-  return lssolve (ComplexMatrix (b), info, rank);
+  double rcond;
+  return lssolve (ComplexMatrix (b), info, rank, rcond);
 }
 
 ComplexMatrix
 ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info) const
 {
   octave_idx_type rank;
-  return lssolve (ComplexMatrix (b), info, rank);
+  double rcond;
+  return lssolve (ComplexMatrix (b), info, rank, rcond);
 }
 
 ComplexMatrix
-ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type& rank) const
+ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info,
+			octave_idx_type& rank) const
 {
-  return lssolve (ComplexMatrix (b), info, rank);
+  double rcond;
+  return lssolve (ComplexMatrix (b), info, rank, rcond);
+}
+
+ComplexMatrix
+ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info,
+			octave_idx_type& rank, double& rcond) const
+{
+  return lssolve (ComplexMatrix (b), info, rank, rcond);
 }
 
 ComplexMatrix
@@ -2416,18 +2427,29 @@
 {
   octave_idx_type info;
   octave_idx_type rank;
-  return lssolve (b, info, rank);
+  double rcond;
+  return lssolve (b, info, rank, rcond);
 }
 
 ComplexMatrix
 ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info) const
 {
   octave_idx_type rank;
-  return lssolve (b, info, rank);
+  double rcond;
+  return lssolve (b, info, rank, rcond);
 }
 
 ComplexMatrix
-ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, octave_idx_type& rank) const
+ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info,
+			octave_idx_type& rank) const
+{
+  double rcond;
+  return lssolve (b, info, rank, rcond);
+}
+
+ComplexMatrix
+ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, 
+			octave_idx_type& rank, double& rcond) const
 {
   ComplexMatrix retval;
 
@@ -2445,7 +2467,7 @@
     {
       volatile octave_idx_type minmn = (m < n ? m : n);
       octave_idx_type maxmn = m > n ? m : n;
-      double rcond = -1.0;
+      rcond = -1.0;
 
       if (m != n)
 	{
@@ -2496,10 +2518,18 @@
 	  if (f77_exception_encountered)
 	    (*current_liboctave_error_handler) 
 	      ("unrecoverable error in zgelsd");
-	  else if (rank < minmn)
-	    (*current_liboctave_warning_handler) 
-	      ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e",
-	       m, n, rank, rcond);
+	  else
+	    {
+	      if (rank < minmn)
+		(*current_liboctave_warning_handler) 
+		  ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e",
+		   m, n, rank, rcond);
+
+	      if (s.elem (0) == 0.0)
+		rcond = 0.0;
+	      else
+		rcond = s.elem (minmn - 1) / s.elem (0);
+	    }
 	}
     }
 
@@ -2511,20 +2541,31 @@
 {
   octave_idx_type info;
   octave_idx_type rank;
-  return lssolve (ComplexColumnVector (b), info, rank);
+  double rcond;
+  return lssolve (ComplexColumnVector (b), info, rank, rcond);
 }
 
 ComplexColumnVector
 ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info) const
 {
   octave_idx_type rank;
-  return lssolve (ComplexColumnVector (b), info, rank);
+  double rcond;
+  return lssolve (ComplexColumnVector (b), info, rank, rcond);
 }
 
 ComplexColumnVector
-ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const
+ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info, 
+			octave_idx_type& rank) const
 {
-  return lssolve (ComplexColumnVector (b), info, rank);
+  double rcond;
+  return lssolve (ComplexColumnVector (b), info, rank, rcond);
+}
+
+ComplexColumnVector
+ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info, 
+			octave_idx_type& rank, double& rcond) const
+{
+  return lssolve (ComplexColumnVector (b), info, rank, rcond);
 }
 
 ComplexColumnVector
@@ -2532,20 +2573,31 @@
 {
   octave_idx_type info;
   octave_idx_type rank;
-  return lssolve (b, info, rank);
+  double rcond;
+  return lssolve (b, info, rank, rcond);
 }
 
 ComplexColumnVector
 ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info) const
 {
   octave_idx_type rank;
-  return lssolve (b, info, rank);
+  double rcond;
+  return lssolve (b, info, rank, rcond);
 }
 
 ComplexColumnVector
 ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info,
 			octave_idx_type& rank) const
 {
+  double rcond;
+  return lssolve (b, info, rank, rcond);
+
+}
+
+ComplexColumnVector
+ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info,
+			octave_idx_type& rank, double& rcond) const
+{
   ComplexColumnVector retval;
 
   octave_idx_type nrhs = 1;
@@ -2562,7 +2614,7 @@
     {
       volatile octave_idx_type minmn = (m < n ? m : n);
       octave_idx_type maxmn = m > n ? m : n;
-      double rcond = -1.0;
+      rcond = -1.0;
 
       if (m != n)
 	{
@@ -2613,9 +2665,17 @@
 	    (*current_liboctave_error_handler) 
 	      ("unrecoverable error in zgelsd");
 	  else if (rank < minmn)
-	    (*current_liboctave_warning_handler) 
-	      ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e",
-	       m, n, rank, rcond);
+	    {
+	      if (rank < minmn)
+		(*current_liboctave_warning_handler) 
+		  ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e",
+		   m, n, rank, rcond);
+
+	      if (s.elem (0) == 0.0)
+		rcond = 0.0;
+	      else
+		rcond = s.elem (minmn - 1) / s.elem (0);
+	    }
 	}
     }