# HG changeset patch
# User jwe
# Date 1193681397 0
# Node ID 0bade2dc44a1d787be84e2fcbbe5350850c79ac7
# Parent  1558d3dab722caacdf8f00292c908dec76276216
[project @ 2007-10-29 18:09:57 by jwe]

diff -r 1558d3dab722 -r 0bade2dc44a1 liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc	Fri Oct 26 18:22:05 2007 +0000
+++ b/liboctave/CMatrix.cc	Mon Oct 29 18:09:57 2007 +0000
@@ -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);
+	    }
 	}
     }
 
diff -r 1558d3dab722 -r 0bade2dc44a1 liboctave/CMatrix.h
--- a/liboctave/CMatrix.h	Fri Oct 26 18:22:05 2007 +0000
+++ b/liboctave/CMatrix.h	Mon Oct 29 18:09:57 2007 +0000
@@ -260,22 +260,35 @@
 
   ComplexMatrix lssolve (const Matrix& b) const;
   ComplexMatrix lssolve (const Matrix& b, octave_idx_type& info) const;
-  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;
+  ComplexMatrix lssolve (const Matrix& b, octave_idx_type& info, 
+			 octave_idx_type& rank, double& rcond) const;
 
   ComplexMatrix lssolve (const ComplexMatrix& b) const;
   ComplexMatrix lssolve (const ComplexMatrix& b, octave_idx_type& info) const;
   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, double& rcond) const;
 
   ComplexColumnVector lssolve (const ColumnVector& b) const;
-  ComplexColumnVector lssolve (const ColumnVector& b, octave_idx_type& info) const;
+  ComplexColumnVector lssolve (const ColumnVector& b,
+			       octave_idx_type& info) const;
   ComplexColumnVector lssolve (const ColumnVector& b, octave_idx_type& info,
 			       octave_idx_type& rank) const;
+  ComplexColumnVector lssolve (const ColumnVector& b, octave_idx_type& info,
+			       octave_idx_type& rank, double& rcond) const;
 
   ComplexColumnVector lssolve (const ComplexColumnVector& b) const;
-  ComplexColumnVector lssolve (const ComplexColumnVector& b, octave_idx_type& info) const;
-  ComplexColumnVector lssolve (const ComplexColumnVector& b, octave_idx_type& info,
+  ComplexColumnVector lssolve (const ComplexColumnVector& b,
+			       octave_idx_type& info) const;
+  ComplexColumnVector lssolve (const ComplexColumnVector& b,
+			       octave_idx_type& info,
 			       octave_idx_type& rank) const;
+  ComplexColumnVector lssolve (const ComplexColumnVector& b,
+			       octave_idx_type& info,
+			       octave_idx_type& rank, double& rcond) const;
 
   ComplexMatrix expm (void) const;
 
diff -r 1558d3dab722 -r 0bade2dc44a1 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Fri Oct 26 18:22:05 2007 +0000
+++ b/liboctave/ChangeLog	Mon Oct 29 18:09:57 2007 +0000
@@ -1,3 +1,42 @@
+2007-10-29  David Bateman  <dbateman@free.fr>
+
+	* CMatrix.h (lssolve (const Matrix&, octave_idx_type&, 
+	octave_idx_type&, double&) const, lssolve (const ComplexMatrix&, 
+	octave_idx_type&, octave_idx_type&, double&) const, lssolve 
+	(const ColumnVector&, octave_idx_type&, octave_idx_type&, 
+	double& rcond) const, lssolve (const ComplexColumnVector&, 
+	octave_idx_type&, octave_idx_type&, double& rcond) const): New
+	declarations.
+	* CMatrix.cc (lssolve (const Matrix&, octave_idx_type&, 
+	octave_idx_type&, double&) const, lssolve (const ComplexMatrix&, 
+	octave_idx_type&, octave_idx_type&, double&) const, lssolve 
+	(const ColumnVector&, octave_idx_type&, octave_idx_type&, 
+	double& rcond) const, lssolve (const ComplexColumnVector&, 
+	octave_idx_type&, octave_idx_type&, double& rcond) const): New
+	methods.
+	(lssolve (const Matrix&, octave_idx_type&, octave_idx_type&,
+	double&) const, lssolve (const ComplexMatrix&, octave_idx_type&, 
+	octave_idx_type&, double&) const): Also return rcond from the
+	singular values returned by XGELSD.
+	* dMatrix.h (lssolve (const Matrix&, octave_idx_type&, 
+	octave_idx_type&, double&) const, lssolve (const ComplexMatrix&, 
+	octave_idx_type&, octave_idx_type&, double&) const, lssolve 
+	(const ColumnVector&, octave_idx_type&, octave_idx_type&, 
+	double& rcond) const, lssolve (const ComplexColumnVector&, 
+	octave_idx_type&, octave_idx_type&, double& rcond) const): New
+	declarations.
+	* dMatrix.cc (lssolve (const Matrix&, octave_idx_type&, 
+	octave_idx_type&, double&) const, lssolve (const ComplexMatrix&, 
+	octave_idx_type&, octave_idx_type&, double&) const, lssolve 
+	(const ColumnVector&, octave_idx_type&, octave_idx_type&, 
+	double& rcond) const, lssolve (const ComplexColumnVector&, 
+	octave_idx_type&, octave_idx_type&, double& rcond) const): New
+	methods.
+	(lssolve (const Matrix&, octave_idx_type&, octave_idx_type&,
+	double&) const, lssolve (const ComplexMatrix&, octave_idx_type&, 
+	octave_idx_type&, double&) const): Also return rcond from the
+	singular values returned by XGELSD.
+		
 2007-10-26  David Bateman  <dbateman@free.fr>
 
 	* dMatrix.cc (Matrix::lssolve): Use xGELSD for rank deficient
diff -r 1558d3dab722 -r 0bade2dc44a1 liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc	Fri Oct 26 18:22:05 2007 +0000
+++ b/liboctave/dMatrix.cc	Mon Oct 29 18:09:57 2007 +0000
@@ -1819,7 +1819,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;
@@ -2039,20 +2039,30 @@
 {
   octave_idx_type info;
   octave_idx_type rank;
-  return lssolve (b, info, rank);
+  double rcond;
+  return lssolve (b, info, rank, rcond);
 }
 
 Matrix
 Matrix::lssolve (const Matrix& b, octave_idx_type& info) const
 {
   octave_idx_type rank;
-  return lssolve (b, info, rank);
+  double rcond;
+  return lssolve (b, info, rank, rcond);
 }
 
 Matrix
 Matrix::lssolve (const Matrix& b, octave_idx_type& info,
 		 octave_idx_type& rank) const
 {
+  double rcond;
+  return lssolve (b, info, rank, rcond);
+}
+
+Matrix
+Matrix::lssolve (const Matrix& b, octave_idx_type& info,
+		 octave_idx_type& rank, double &rcond) const
+{
   Matrix retval;
 
   octave_idx_type nrhs = b.cols ();
@@ -2069,7 +2079,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)
 	{
 	  retval = Matrix (maxmn, nrhs, 0.0);
@@ -2118,9 +2128,16 @@
 	  if (f77_exception_encountered)
 	    (*current_liboctave_error_handler) 
 	      ("unrecoverable error in dgelsd");
-	  else if (rank < minmn)
-	    (*current_liboctave_warning_handler) 
-	      ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank);
+	  else 
+	    {
+	      if (rank < minmn)
+		(*current_liboctave_warning_handler) 
+		  ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank);
+	      if (s.elem (0) == 0.0)
+		rcond = 0.0;
+	      else
+		rcond = s.elem (minmn - 1) / s.elem (0);
+	    }
 	}
     }
 
@@ -2133,7 +2150,8 @@
   ComplexMatrix tmp (*this);
   octave_idx_type info;
   octave_idx_type rank;
-  return tmp.lssolve (b, info, rank);
+  double rcond;
+  return tmp.lssolve (b, info, rank, rcond);
 }
 
 ComplexMatrix
@@ -2141,14 +2159,25 @@
 {
   ComplexMatrix tmp (*this);
   octave_idx_type rank;
-  return tmp.lssolve (b, info, rank);
+  double rcond;
+  return tmp.lssolve (b, info, rank, rcond);
 }
 
 ComplexMatrix
-Matrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, octave_idx_type& rank) const
+Matrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, 
+		 octave_idx_type& rank) const
 {
   ComplexMatrix tmp (*this);
-  return tmp.lssolve (b, info, rank);
+  double rcond;
+  return tmp.lssolve (b, info, rank, rcond);
+}
+
+ComplexMatrix
+Matrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, 
+		 octave_idx_type& rank, double& rcond) const
+{
+  ComplexMatrix tmp (*this);
+  return tmp.lssolve (b, info, rank, rcond);
 }
 
 ColumnVector
@@ -2156,20 +2185,30 @@
 {
   octave_idx_type info;
   octave_idx_type rank;
-  return lssolve (b, info, rank);
+  double rcond;
+  return lssolve (b, info, rank, rcond);
 }
 
 ColumnVector
 Matrix::lssolve (const ColumnVector& b, octave_idx_type& info) const
 {
   octave_idx_type rank;
-  return lssolve (b, info, rank);
+  double rcond;
+  return lssolve (b, info, rank, rcond);
 }
 
 ColumnVector
 Matrix::lssolve (const ColumnVector& b, octave_idx_type& info,
 		 octave_idx_type& rank) const
 {
+  double rcond;
+  return lssolve (b, info, rank, rcond);
+}
+
+ColumnVector
+Matrix::lssolve (const ColumnVector& b, octave_idx_type& info,
+		 octave_idx_type& rank, double &rcond) const
+{
   ColumnVector retval;
 
   octave_idx_type nrhs = 1;
@@ -2186,7 +2225,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)
 	{
@@ -2236,8 +2275,15 @@
 	    (*current_liboctave_error_handler) 
 	      ("unrecoverable error in dgelsd");
 	  else if (rank < minmn)
-	    (*current_liboctave_warning_handler) 
-	      ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank);
+	    {
+	      if (rank < minmn)
+		(*current_liboctave_warning_handler) 
+		  ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank);
+	      if (s.elem (0) == 0.0)
+		rcond = 0.0;
+	      else
+		rcond = s.elem (minmn - 1) / s.elem (0);
+	    }
 	}
     }
 
@@ -2248,21 +2294,36 @@
 Matrix::lssolve (const ComplexColumnVector& b) const
 {
   ComplexMatrix tmp (*this);
-  return tmp.lssolve (b);
+  octave_idx_type info;
+  octave_idx_type rank;
+  double rcond;
+  return tmp.lssolve (b, info, rank, rcond);
 }
 
 ComplexColumnVector
 Matrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info) const
 {
   ComplexMatrix tmp (*this);
-  return tmp.lssolve (b, info);
+  octave_idx_type rank;
+  double rcond;
+  return tmp.lssolve (b, info, rank, rcond);
 }
 
 ComplexColumnVector
-Matrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const
+Matrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, 
+		 octave_idx_type& rank) const
 {
   ComplexMatrix tmp (*this);
-  return tmp.lssolve (b, info, rank);
+  double rcond;
+  return tmp.lssolve (b, info, rank, rcond);
+}
+
+ComplexColumnVector
+Matrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, 
+		 octave_idx_type& rank, double &rcond) const
+{
+  ComplexMatrix tmp (*this);
+  return tmp.lssolve (b, info, rank, rcond);
 }
 
 // Constants for matrix exponential calculation.
diff -r 1558d3dab722 -r 0bade2dc44a1 liboctave/dMatrix.h
--- a/liboctave/dMatrix.h	Fri Oct 26 18:22:05 2007 +0000
+++ b/liboctave/dMatrix.h	Mon Oct 29 18:09:57 2007 +0000
@@ -226,21 +226,34 @@
   // Singular solvers
   Matrix lssolve (const Matrix& b) const;
   Matrix lssolve (const Matrix& b, octave_idx_type& info) const;
-  Matrix lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type& rank) const;
+  Matrix lssolve (const Matrix& b, octave_idx_type& info, 
+		  octave_idx_type& rank) const;
+  Matrix lssolve (const Matrix& b, octave_idx_type& info, 
+		  octave_idx_type& rank, double& rcond) const;
 
   ComplexMatrix lssolve (const ComplexMatrix& b) const;
   ComplexMatrix lssolve (const ComplexMatrix& b, octave_idx_type& info) const;
   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, double &rcond) const;
 
   ColumnVector lssolve (const ColumnVector& b) const;
   ColumnVector lssolve (const ColumnVector& b, octave_idx_type& info) const;
-  ColumnVector lssolve (const ColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const;
+  ColumnVector lssolve (const ColumnVector& b, octave_idx_type& info,
+			octave_idx_type& rank) const;
+  ColumnVector lssolve (const ColumnVector& b, octave_idx_type& info,
+			octave_idx_type& rank, double& rcond) const;
 
   ComplexColumnVector lssolve (const ComplexColumnVector& b) const;
-  ComplexColumnVector lssolve (const ComplexColumnVector& b, octave_idx_type& info) const;
-  ComplexColumnVector lssolve (const ComplexColumnVector& b, octave_idx_type& info,
+  ComplexColumnVector lssolve (const ComplexColumnVector& b, 
+			       octave_idx_type& info) const;
+  ComplexColumnVector lssolve (const ComplexColumnVector& b,
+			       octave_idx_type& info,
 			       octave_idx_type& rank) const;
+  ComplexColumnVector lssolve (const ComplexColumnVector& b, 
+			       octave_idx_type& info,
+			       octave_idx_type& rank, double& rcond) const;
 
   Matrix expm (void) const;