Mercurial > hg > octave-terminal
changeset 7076:0bade2dc44a1
[project @ 2007-10-29 18:09:57 by jwe]
author | jwe |
---|---|
date | Mon, 29 Oct 2007 18:09:57 +0000 |
parents | 1558d3dab722 |
children | 525cd5f47ab6 |
files | liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/ChangeLog liboctave/dMatrix.cc liboctave/dMatrix.h |
diffstat | 5 files changed, 237 insertions(+), 51 deletions(-) [+] |
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); + } } }
--- a/liboctave/CMatrix.h +++ b/liboctave/CMatrix.h @@ -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;
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -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
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -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.
--- a/liboctave/dMatrix.h +++ b/liboctave/dMatrix.h @@ -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;