# HG changeset patch # User dbateman # Date 1142956908 0 # Node ID 2fe20065a545777c62b9e80455bac8f9b47a82ea # Parent 70cc04f9af41f71434919ab6cbf3ba2e917452b8 [project @ 2006-03-21 16:01:46 by dbateman] diff --git a/liboctave/CSparse.cc b/liboctave/CSparse.cc --- a/liboctave/CSparse.cc +++ b/liboctave/CSparse.cc @@ -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 diff --git a/liboctave/CSparse.h b/liboctave/CSparse.h --- a/liboctave/CSparse.h +++ b/liboctave/CSparse.h @@ -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& vec, bool inv = false) const; + SparseComplexMatrix permute (const Array& vec, + bool inv = false) const; SparseComplexMatrix ipermute (const Array& vec) const; diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog --- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,14 @@ +2006-03-21 David Bateman + + * 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 * str-vec.cc (vector::list_in_columns): New optional arg, width. diff --git a/liboctave/dSparse.cc b/liboctave/dSparse.cc --- a/liboctave/dSparse.cc +++ b/liboctave/dSparse.cc @@ -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 diff --git a/liboctave/dSparse.h b/liboctave/dSparse.h --- a/liboctave/dSparse.h +++ b/liboctave/dSparse.h @@ -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, diff --git a/liboctave/sparse-dmsolve.cc b/liboctave/sparse-dmsolve.cc --- a/liboctave/sparse-dmsolve.cc +++ b/liboctave/sparse-dmsolve.cc @@ -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;