changeset 5697:2fe20065a545

[project @ 2006-03-21 16:01:46 by dbateman]
author dbateman
date Tue, 21 Mar 2006 16:01:48 +0000
parents 70cc04f9af41
children e33aff8ba378
files liboctave/CSparse.cc liboctave/CSparse.h liboctave/ChangeLog liboctave/dSparse.cc liboctave/dSparse.h liboctave/sparse-dmsolve.cc
diffstat 6 files changed, 117 insertions(+), 76 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CSparse.cc	Mon Mar 20 18:52:46 2006 +0000
+++ b/liboctave/CSparse.cc	Tue Mar 21 16:01:48 2006 +0000
@@ -6717,16 +6717,17 @@
 }
 
 ComplexMatrix
-SparseComplexMatrix::solve (SparseType &mattype, const Matrix& b, octave_idx_type& info, 
-			    double& rcond) const
+SparseComplexMatrix::solve (SparseType &mattype, const Matrix& b,
+			    octave_idx_type& info, double& rcond) const
 {
   return solve (mattype, b, info, rcond, 0);
 }
 
 ComplexMatrix
-SparseComplexMatrix::solve (SparseType &mattype, const Matrix& b, octave_idx_type& err, 
-			    double& rcond, 
-			    solve_singularity_handler sing_handler) const
+SparseComplexMatrix::solve (SparseType &mattype, const Matrix& b,
+			    octave_idx_type& err, double& rcond, 
+			    solve_singularity_handler sing_handler,
+			    bool singular_fallback) const
 {
   ComplexMatrix retval;
   int typ = mattype.type (false);
@@ -6753,7 +6754,7 @@
       return ComplexMatrix ();
     }
 
-  if (mattype.type(false) == SparseType::Rectangular)
+  if (singular_fallback && mattype.type(false) == SparseType::Rectangular)
     {
       rcond = 1.;
 #ifdef USE_QRSOLVE
@@ -6793,7 +6794,8 @@
 SparseComplexMatrix
 SparseComplexMatrix::solve (SparseType &mattype, const SparseMatrix& b, 
 			    octave_idx_type& err, double& rcond,
-			    solve_singularity_handler sing_handler) const
+			    solve_singularity_handler sing_handler,
+			    bool singular_fallback) const
 {
   SparseComplexMatrix retval;
   int typ = mattype.type (false);
@@ -6820,7 +6822,7 @@
       return SparseComplexMatrix ();
     }
 
-  if (mattype.type(false) == SparseType::Rectangular)
+  if (singular_fallback && mattype.type(false) == SparseType::Rectangular)
     {
       rcond = 1.;
 #ifdef USE_QRSOLVE
@@ -6852,15 +6854,16 @@
 
 ComplexMatrix
 SparseComplexMatrix::solve (SparseType &mattype, const ComplexMatrix& b, 
-		     octave_idx_type& info, double& rcond) const
+			    octave_idx_type& info, double& rcond) const
 {
   return solve (mattype, b, info, rcond, 0);
 }
 
 ComplexMatrix
 SparseComplexMatrix::solve (SparseType &mattype, const ComplexMatrix& b, 
-		     octave_idx_type& err, double& rcond, 
-		     solve_singularity_handler sing_handler) const
+			    octave_idx_type& err, double& rcond, 
+			    solve_singularity_handler sing_handler,
+			    bool singular_fallback) const
 {
   ComplexMatrix retval;
   int typ = mattype.type (false);
@@ -6887,7 +6890,7 @@
       return ComplexMatrix ();
     }
 
-  if (mattype.type(false) == SparseType::Rectangular)
+  if (singular_fallback && mattype.type(false) == SparseType::Rectangular)
     {
       rcond = 1.;
 #ifdef USE_QRSOLVE
@@ -6912,7 +6915,7 @@
 
 SparseComplexMatrix
 SparseComplexMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b, 
-		     octave_idx_type& info) const
+			    octave_idx_type& info) const
 {
   double rcond;
   return solve (mattype, b, info, rcond, 0);
@@ -6920,7 +6923,7 @@
 
 SparseComplexMatrix
 SparseComplexMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b,
-		     octave_idx_type& info, double& rcond) const
+			    octave_idx_type& info, double& rcond) const
 {
   return solve (mattype, b, info, rcond, 0);
 }
@@ -6928,7 +6931,8 @@
 SparseComplexMatrix
 SparseComplexMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b, 
 			    octave_idx_type& err, double& rcond,
-			    solve_singularity_handler sing_handler) const
+			    solve_singularity_handler sing_handler,
+			    bool singular_fallback) const
 {
   SparseComplexMatrix retval;
   int typ = mattype.type (false);
@@ -6955,7 +6959,7 @@
       return SparseComplexMatrix ();
     }
 
-  if (mattype.type(false) == SparseType::Rectangular)
+  if (singular_fallback && mattype.type(false) == SparseType::Rectangular)
     {
       rcond = 1.;
 #ifdef USE_QRSOLVE
--- a/liboctave/CSparse.h	Mon Mar 20 18:52:46 2006 +0000
+++ b/liboctave/CSparse.h	Tue Mar 21 16:01:48 2006 +0000
@@ -273,27 +273,33 @@
 public:
   // Generic interface to solver with no probing of type
   ComplexMatrix solve (SparseType &typ, const Matrix& b) const;
-  ComplexMatrix solve (SparseType &typ, const Matrix& b, octave_idx_type& info) const;
+  ComplexMatrix solve (SparseType &typ, const Matrix& b, 
+		       octave_idx_type& info) const;
   ComplexMatrix solve (SparseType &typ, const Matrix& b, octave_idx_type& info, 
-		double& rcond) const;
+		       double& rcond) const;
   ComplexMatrix solve (SparseType &typ, const Matrix& b, octave_idx_type& info, 
-		double& rcond, solve_singularity_handler sing_handler) const;
+		       double& rcond, solve_singularity_handler sing_handler,
+		       bool singular_fallback = true) const;
 
   ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b) const;
   ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b, 
 		       octave_idx_type& info) const;
-  ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b, octave_idx_type& info, 
-		       double& rcond) const;
-  ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b, octave_idx_type& info, 
-		double& rcond, solve_singularity_handler sing_handler) const;
+  ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b, 
+		       octave_idx_type& info, double& rcond) const;
+  ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b,
+		       octave_idx_type& info, double& rcond, 
+		       solve_singularity_handler sing_handler,
+		       bool singular_fallback = true) const;
 
   SparseComplexMatrix solve (SparseType &typ, const SparseMatrix& b) const;
   SparseComplexMatrix solve (SparseType &typ, const SparseMatrix& b, 
-		      octave_idx_type& info) const;
-  SparseComplexMatrix solve (SparseType &typ, const SparseMatrix& b, octave_idx_type& info, 
-		      double& rcond) const;
-  SparseComplexMatrix solve (SparseType &typ, const SparseMatrix& b, octave_idx_type& info, 
-		double& rcond, solve_singularity_handler sing_handler) const;
+			     octave_idx_type& info) const;
+  SparseComplexMatrix solve (SparseType &typ, const SparseMatrix& b, 
+			     octave_idx_type& info, double& rcond) const;
+  SparseComplexMatrix solve (SparseType &typ, const SparseMatrix& b, 
+			     octave_idx_type& info, double& rcond, 
+			     solve_singularity_handler sing_handler,
+			     bool singular_fallback = true) const;
 
   SparseComplexMatrix solve (SparseType &typ, 
 			     const SparseComplexMatrix& b) const;
@@ -301,21 +307,24 @@
 			     octave_idx_type& info) const;
   SparseComplexMatrix solve (SparseType &typ, const SparseComplexMatrix& b, 
 			     octave_idx_type& info, double& rcond) const;
-  SparseComplexMatrix solve (SparseType &typ, const SparseComplexMatrix& b, octave_idx_type& info, 
-	       double& rcond, solve_singularity_handler sing_handler) const;
+  SparseComplexMatrix solve (SparseType &typ, const SparseComplexMatrix& b,
+			     octave_idx_type& info, double& rcond, 
+			     solve_singularity_handler sing_handler,
+			     bool singular_fallback = true) const;
 
   ComplexColumnVector solve (SparseType &typ, const ColumnVector& b) const;
   ComplexColumnVector solve (SparseType &typ, const ColumnVector& b, 
-		      octave_idx_type& info) const;
+			     octave_idx_type& info) const;
   ComplexColumnVector solve (SparseType &typ, const ColumnVector& b, 
-		      octave_idx_type& info, double& rcond) const;
-  ComplexColumnVector solve (SparseType &typ, const ColumnVector& b, octave_idx_type& info,
-		double& rcond, solve_singularity_handler sing_handler) const;
+			     octave_idx_type& info, double& rcond) const;
+  ComplexColumnVector solve (SparseType &typ, const ColumnVector& b,
+			     octave_idx_type& info, double& rcond, 
+			     solve_singularity_handler sing_handler) const;
 
   ComplexColumnVector solve (SparseType &typ, 
 			     const ComplexColumnVector& b) const;
-  ComplexColumnVector solve (SparseType &typ, 
-			     const ComplexColumnVector& b, octave_idx_type& info) const;
+  ComplexColumnVector solve (SparseType &typ, const ComplexColumnVector& b, 
+			     octave_idx_type& info) const;
   ComplexColumnVector solve (SparseType &typ, const ComplexColumnVector& b,
 			     octave_idx_type& info, double& rcond) const;
   ComplexColumnVector solve (SparseType &typ, const ComplexColumnVector& b,
@@ -325,7 +334,8 @@
   // Generic interface to solver with probing of type
   ComplexMatrix solve (const Matrix& b) const;
   ComplexMatrix solve (const Matrix& b, octave_idx_type& info) const;
-  ComplexMatrix solve (const Matrix& b, octave_idx_type& info, double& rcond) const;
+  ComplexMatrix solve (const Matrix& b, octave_idx_type& info, 
+		       double& rcond) const;
   ComplexMatrix solve (const Matrix& b, octave_idx_type& info, double& rcond, 
 		       solve_singularity_handler sing_handler) const;
 
@@ -333,8 +343,8 @@
   ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info) const;
   ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, 
 		       double& rcond) const;
-  ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, double& rcond,
-		       solve_singularity_handler sing_handler) const;
+  ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info,
+		       double& rcond, solve_singularity_handler sing_handler) const;
 
   SparseComplexMatrix solve (const SparseMatrix& b) const;
   SparseComplexMatrix solve (const SparseMatrix& b, octave_idx_type& info) const;
@@ -342,25 +352,28 @@
 			     double& rcond) const;
   SparseComplexMatrix solve (const SparseMatrix& b, octave_idx_type& info, 
 			     double& rcond, 
-		       solve_singularity_handler sing_handler) const;
+			     solve_singularity_handler sing_handler) const;
 
   SparseComplexMatrix solve (const SparseComplexMatrix& b) const;
-  SparseComplexMatrix solve (const SparseComplexMatrix& b, octave_idx_type& info) const;
-  SparseComplexMatrix solve (const SparseComplexMatrix& b, octave_idx_type& info, 
-			     double& rcond) const;
-  SparseComplexMatrix solve (const SparseComplexMatrix& b, octave_idx_type& info, 
-			     double& rcond,
+  SparseComplexMatrix solve (const SparseComplexMatrix& b, 
+			     octave_idx_type& info) const;
+  SparseComplexMatrix solve (const SparseComplexMatrix& b,
+			     octave_idx_type& info, double& rcond) const;
+  SparseComplexMatrix solve (const SparseComplexMatrix& b,
+			     octave_idx_type& info, double& rcond,
 			     solve_singularity_handler sing_handler) const;
 
   ComplexColumnVector solve (const ColumnVector& b) const;
   ComplexColumnVector solve (const ColumnVector& b, octave_idx_type& info) const;
   ComplexColumnVector solve (const ColumnVector& b, octave_idx_type& info,
 			     double& rcond) const;
-  ComplexColumnVector solve (const ColumnVector& b, octave_idx_type& info, double& rcond,
+  ComplexColumnVector solve (const ColumnVector& b, octave_idx_type& info, 
+			     double& rcond,
 			     solve_singularity_handler sing_handler) const;
 
   ComplexColumnVector solve (const ComplexColumnVector& b) const;
-  ComplexColumnVector solve (const ComplexColumnVector& b, octave_idx_type& info) const;
+  ComplexColumnVector solve (const ComplexColumnVector& b, 
+			     octave_idx_type& info) const;
   ComplexColumnVector solve (const ComplexColumnVector& b, octave_idx_type& info,
 			     double& rcond) const;
   ComplexColumnVector solve (const ComplexColumnVector& b, octave_idx_type& info,
@@ -377,7 +390,8 @@
 
   SparseComplexMatrix reshape (const dim_vector& new_dims) const;
 
-  SparseComplexMatrix permute (const Array<octave_idx_type>& vec, bool inv = false) const;
+  SparseComplexMatrix permute (const Array<octave_idx_type>& vec, 
+			       bool inv = false) const;
 
   SparseComplexMatrix ipermute (const Array<octave_idx_type>& vec) const;
 
--- a/liboctave/ChangeLog	Mon Mar 20 18:52:46 2006 +0000
+++ b/liboctave/ChangeLog	Tue Mar 21 16:01:48 2006 +0000
@@ -1,3 +1,14 @@
+2006-03-21  David Bateman  <dbateman@free.fr>
+
+	    * dSparse.cc (solve): Add argument singular_fallback, to allow
+	    fallback to QR solvers to be optional.
+	    * CSparse.cc (solve): ditto.
+	    * dSparse.h (solve): update declaration for new argument.
+	    * CSparse.h (solve): ditto.
+	    * sparse-dmsolve.cc (dmsolve): Use singular_fallback argument
+	    to bypass QR solvers when solving the well determined part of
+	    the problem.
+
 2006-03-17  John W. Eaton  <jwe@octave.org>
 
 	* str-vec.cc (vector::list_in_columns): New optional arg, width.
--- a/liboctave/dSparse.cc	Mon Mar 20 18:52:46 2006 +0000
+++ b/liboctave/dSparse.cc	Tue Mar 21 16:01:48 2006 +0000
@@ -6884,7 +6884,8 @@
 }
 
 Matrix
-SparseMatrix::solve (SparseType &mattype, const Matrix& b, octave_idx_type& info) const
+SparseMatrix::solve (SparseType &mattype, const Matrix& b, 
+		     octave_idx_type& info) const
 {
   double rcond;
   return solve (mattype, b, info, rcond, 0);
@@ -6899,8 +6900,8 @@
 
 Matrix
 SparseMatrix::solve (SparseType &mattype, const Matrix& b, octave_idx_type& err, 
-		     double& rcond, 
-		     solve_singularity_handler sing_handler) const
+		     double& rcond, solve_singularity_handler sing_handler,
+		     bool singular_fallback) const
 {
   Matrix retval;
   int typ = mattype.type (false);
@@ -6929,7 +6930,7 @@
     }
 
   // Rectangular or one of the above solvers flags a singular matrix
-  if (mattype.type (false) == SparseType::Rectangular)
+  if (singular_fallback && mattype.type (false) == SparseType::Rectangular)
     {
       rcond = 1.;
 #ifdef USE_QRSOLVE
@@ -6968,7 +6969,8 @@
 SparseMatrix
 SparseMatrix::solve (SparseType &mattype, const SparseMatrix& b, 
 		     octave_idx_type& err, double& rcond,
-		     solve_singularity_handler sing_handler) const
+		     solve_singularity_handler sing_handler,
+		     bool singular_fallback) const
 {
   SparseMatrix retval;
   int typ = mattype.type (false);
@@ -6995,7 +6997,7 @@
       return SparseMatrix ();
     }
 
-  if (mattype.type (false) == SparseType::Rectangular)
+  if (singular_fallback && mattype.type (false) == SparseType::Rectangular)
     {
       rcond = 1.;
 #ifdef USE_QRSOLVE
@@ -7035,7 +7037,8 @@
 ComplexMatrix
 SparseMatrix::solve (SparseType &mattype, const ComplexMatrix& b, 
 		     octave_idx_type& err, double& rcond, 
-		     solve_singularity_handler sing_handler) const
+		     solve_singularity_handler sing_handler,
+		     bool singular_fallback) const
 {
   ComplexMatrix retval;
   int typ = mattype.type (false);
@@ -7062,7 +7065,7 @@
       return ComplexMatrix ();
     }
 
-  if (mattype.type(false) == SparseType::Rectangular)
+  if (singular_fallback && mattype.type(false) == SparseType::Rectangular)
     {
       rcond = 1.;
 #ifdef USE_QRSOLVE
@@ -7102,7 +7105,8 @@
 SparseComplexMatrix
 SparseMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b, 
 		     octave_idx_type& err, double& rcond,
-		     solve_singularity_handler sing_handler) const
+		     solve_singularity_handler sing_handler,
+		     bool singular_fallback) const
 {
   SparseComplexMatrix retval;
   int typ = mattype.type (false);
@@ -7129,7 +7133,7 @@
       return SparseComplexMatrix ();
     }
 
-  if (mattype.type(false) == SparseType::Rectangular)
+  if (singular_fallback && mattype.type(false) == SparseType::Rectangular)
     {
       rcond = 1.;
 #ifdef USE_QRSOLVE
--- a/liboctave/dSparse.h	Mon Mar 20 18:52:46 2006 +0000
+++ b/liboctave/dSparse.h	Tue Mar 21 16:01:48 2006 +0000
@@ -265,24 +265,29 @@
   Matrix solve (SparseType &typ, const Matrix& b, octave_idx_type& info) const;
   Matrix solve (SparseType &typ, const Matrix& b, octave_idx_type& info, 
 		double& rcond) const;
-  Matrix solve (SparseType &typ, const Matrix& b, octave_idx_type& info, double& rcond,
-		solve_singularity_handler sing_handler) const;
+  Matrix solve (SparseType &typ, const Matrix& b, octave_idx_type& info,
+		double& rcond, solve_singularity_handler sing_handler,
+		bool singular_fallback = true) const;
 
   ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b) const;
   ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b, 
 		       octave_idx_type& info) const;
-  ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b, octave_idx_type& info, 
-		       double& rcond) const;
-  ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b, octave_idx_type& info, 
-		double& rcond, solve_singularity_handler sing_handler) const;
+  ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b, 
+		       octave_idx_type& info, double& rcond) const;
+  ComplexMatrix solve (SparseType &typ, const ComplexMatrix& b,
+		       octave_idx_type& info, double& rcond, 
+		       solve_singularity_handler sing_handler,
+		       bool singular_fallback = true) const;
 
   SparseMatrix solve (SparseType &typ, const SparseMatrix& b) const;
   SparseMatrix solve (SparseType &typ, const SparseMatrix& b, 
 		      octave_idx_type& info) const;
-  SparseMatrix solve (SparseType &typ, const SparseMatrix& b, octave_idx_type& info, 
-		      double& rcond) const;
-  SparseMatrix solve (SparseType &typ, const SparseMatrix& b, octave_idx_type& info, 
-		double& rcond, solve_singularity_handler sing_handler) const;
+  SparseMatrix solve (SparseType &typ, const SparseMatrix& b,
+		      octave_idx_type& info, double& rcond) const;
+  SparseMatrix solve (SparseType &typ, const SparseMatrix& b,
+		      octave_idx_type& info, double& rcond, 
+		      solve_singularity_handler sing_handler,
+		      bool singular_fallback = true) const;
 
   SparseComplexMatrix solve (SparseType &typ, 
 			     const SparseComplexMatrix& b) const;
@@ -290,21 +295,24 @@
 			     octave_idx_type& info) const;
   SparseComplexMatrix solve (SparseType &typ, const SparseComplexMatrix& b, 
 			     octave_idx_type& info, double& rcond) const;
-  SparseComplexMatrix solve (SparseType &typ, const SparseComplexMatrix& b, octave_idx_type& info, 
-	       double& rcond, solve_singularity_handler sing_handler) const;
+  SparseComplexMatrix solve (SparseType &typ, const SparseComplexMatrix& b,
+			     octave_idx_type& info, double& rcond, 
+			     solve_singularity_handler sing_handler,
+			     bool singular_fallabck = true) const;
 
   ColumnVector solve (SparseType &typ, const ColumnVector& b) const;
   ColumnVector solve (SparseType &typ, const ColumnVector& b, 
 		      octave_idx_type& info) const;
   ColumnVector solve (SparseType &typ, const ColumnVector& b, 
 		      octave_idx_type& info, double& rcond) const;
-  ColumnVector solve (SparseType &typ, const ColumnVector& b, octave_idx_type& info,
-		double& rcond, solve_singularity_handler sing_handler) const;
+  ColumnVector solve (SparseType &typ, const ColumnVector& b,
+		      octave_idx_type& info, double& rcond, 
+		      solve_singularity_handler sing_handler) const;
 
   ComplexColumnVector solve (SparseType &typ, 
 			     const ComplexColumnVector& b) const;
-  ComplexColumnVector solve (SparseType &typ, 
-			     const ComplexColumnVector& b, octave_idx_type& info) const;
+  ComplexColumnVector solve (SparseType &typ, const ComplexColumnVector& b, 
+			     octave_idx_type& info) const;
   ComplexColumnVector solve (SparseType &typ, const ComplexColumnVector& b,
 			     octave_idx_type& info, double& rcond) const;
   ComplexColumnVector solve (SparseType &typ, const ComplexColumnVector& b,
--- a/liboctave/sparse-dmsolve.cc	Mon Mar 20 18:52:46 2006 +0000
+++ b/liboctave/sparse-dmsolve.cc	Tue Mar 21 16:01:48 2006 +0000
@@ -395,7 +395,6 @@
 	pinv [p [i]] = i;
       RT btmp;
       dmsolve_permute (btmp, b, pinv);
-      SparseType mtyp (SparseType::Full);
       info = 0;
       retval.resize (nc, b_nc);
 
@@ -431,8 +430,9 @@
 	  RT btmp2 = dmsolve_extract (btmp, NULL, NULL, dm->rr [1], dm->rr [2], 
 				      0, b_nc);
 	  double rcond = 0.0;
+	  SparseType mtyp (SparseType::Full);
 	  RT mtmp = m.solve (mtyp, btmp2, info, rcond, 
-			     solve_singularity_warning);	
+			     solve_singularity_warning, false);	
 	  if (info != 0)
 	    {
 	      info = 0;