Mercurial > octave-nkf
changeset 7788:45f5faba05a2
Add the rcond function
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Wed, 14 May 2008 18:09:56 +0200 |
| parents | 6b521b1e3631 |
| children | 82be108cc558 |
| files | liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/ChangeLog liboctave/dMatrix.cc liboctave/dMatrix.h src/ChangeLog src/DLD-FUNCTIONS/rcond.cc src/Makefile.in |
| diffstat | 8 files changed, 821 insertions(+), 389 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc Tue May 13 21:12:12 2008 +0200 +++ b/liboctave/CMatrix.cc Wed May 14 18:09:56 2008 +0200 @@ -983,45 +983,45 @@ ComplexMatrix::inverse (void) const { octave_idx_type info; - double rcond; + double rcon; MatrixType mattype (*this); - return inverse (mattype, info, rcond, 0, 0); + return inverse (mattype, info, rcon, 0, 0); } ComplexMatrix ComplexMatrix::inverse (octave_idx_type& info) const { - double rcond; + double rcon; MatrixType mattype (*this); - return inverse (mattype, info, rcond, 0, 0); + return inverse (mattype, info, rcon, 0, 0); } ComplexMatrix -ComplexMatrix::inverse (octave_idx_type& info, double& rcond, int force, +ComplexMatrix::inverse (octave_idx_type& info, double& rcon, int force, int calc_cond) const { MatrixType mattype (*this); - return inverse (mattype, info, rcond, force, calc_cond); + return inverse (mattype, info, rcon, force, calc_cond); } ComplexMatrix ComplexMatrix::inverse (MatrixType &mattype) const { octave_idx_type info; - double rcond; - return inverse (mattype, info, rcond, 0, 0); + double rcon; + return inverse (mattype, info, rcon, 0, 0); } ComplexMatrix ComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info) const { - double rcond; - return inverse (mattype, info, rcond, 0, 0); + double rcon; + return inverse (mattype, info, rcon, 0, 0); } ComplexMatrix ComplexMatrix::tinverse (MatrixType &mattype, octave_idx_type& info, - double& rcond, int force, int calc_cond) const + double& rcon, int force, int calc_cond) const { ComplexMatrix retval; @@ -1045,7 +1045,7 @@ F77_CHAR_ARG_LEN (1))); // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + rcon = 0.0; if (info != 0) info = -1; else if (calc_cond) @@ -1059,7 +1059,7 @@ F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&job, 1), F77_CONST_CHAR_ARG2 (&uplo, 1), F77_CONST_CHAR_ARG2 (&udiag, 1), - nr, tmp_data, nr, rcond, + nr, tmp_data, nr, rcon, cwork, rwork, ztrcon_info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) @@ -1078,7 +1078,7 @@ ComplexMatrix ComplexMatrix::finverse (MatrixType &mattype, octave_idx_type& info, - double& rcond, int force, int calc_cond) const + double& rcon, int force, int calc_cond) const { ComplexMatrix retval; @@ -1118,7 +1118,7 @@ F77_XFCN (zgetrf, ZGETRF, (nc, nc, tmp_data, nr, pipvt, info)); // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + rcon = 0.0; if (info != 0) info = -1; else if (calc_cond) @@ -1130,7 +1130,7 @@ double *prz = rz.fortran_vec (); F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, tmp_data, nr, anorm, - rcond, pz, prz, zgecon_info + rcon, pz, prz, zgecon_info F77_CHAR_ARG_LEN (1))); if (zgecon_info != 0) @@ -1159,7 +1159,7 @@ ComplexMatrix ComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info, - double& rcond, int force, int calc_cond) const + double& rcon, int force, int calc_cond) const { int typ = mattype.type (false); ComplexMatrix ret; @@ -1168,7 +1168,7 @@ typ = mattype.type (*this); if (typ == MatrixType::Upper || typ == MatrixType::Lower) - ret = tinverse (mattype, info, rcond, force, calc_cond); + ret = tinverse (mattype, info, rcon, force, calc_cond); else { if (mattype.is_hermitian ()) @@ -1177,9 +1177,9 @@ if (info == 0) { if (calc_cond) - rcond = chol.rcond(); + rcon = chol.rcond(); else - rcond = 1.0; + rcon = 1.0; ret = chol.inverse (); } else @@ -1187,9 +1187,9 @@ } if (!mattype.is_hermitian ()) - ret = finverse(mattype, info, rcond, force, calc_cond); - - if ((mattype.is_hermitian () || calc_cond) && rcond == 0.) + ret = finverse(mattype, info, rcon, force, calc_cond); + + if ((mattype.is_hermitian () || calc_cond) && rcon == 0.) ret = ComplexMatrix (rows (), columns (), Complex (octave_Inf, 0.)); } @@ -1553,19 +1553,19 @@ ComplexMatrix::determinant (void) const { octave_idx_type info; - double rcond; - return determinant (info, rcond, 0); + double rcon; + return determinant (info, rcon, 0); } ComplexDET ComplexMatrix::determinant (octave_idx_type& info) const { - double rcond; - return determinant (info, rcond, 0); + double rcon; + return determinant (info, rcon, 0); } ComplexDET -ComplexMatrix::determinant (octave_idx_type& info, double& rcond, int calc_cond) const +ComplexMatrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const { ComplexDET retval; @@ -1594,7 +1594,7 @@ F77_XFCN (zgetrf, ZGETRF, (nr, nc, tmp_data, nr, pipvt, info)); // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + rcon = 0.0; if (info != 0) { info = -1; @@ -1613,7 +1613,7 @@ F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, tmp_data, nr, anorm, - rcond, pz, prz, info + rcon, pz, prz, info F77_CHAR_ARG_LEN (1))); } @@ -1658,9 +1658,176 @@ return retval; } +double +ComplexMatrix::rcond (void) const +{ + MatrixType mattype (*this); + return rcond (mattype); +} + +double +ComplexMatrix::rcond (MatrixType &mattype) const +{ + double rcon; + octave_idx_type nr = rows (); + octave_idx_type nc = cols (); + + if (nr != nc) + (*current_liboctave_error_handler) ("matrix must be square"); + else if (nr == 0 || nc == 0) + rcon = octave_Inf; + else + { + int typ = mattype.type (); + + if (typ == MatrixType::Unknown) + typ = mattype.type (*this); + + // Only calculate the condition number for LU/Cholesky + if (typ == MatrixType::Upper) + { + const Complex *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'U'; + char dia = 'N'; + + Array<Complex> z (2 * nc); + Complex *pz = z.fortran_vec (); + Array<double> rz (nc); + double *prz = rz.fortran_vec (); + + F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, prz, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0; + } + else if (typ == MatrixType::Permuted_Upper) + (*current_liboctave_error_handler) + ("permuted triangular matrix not implemented"); + else if (typ == MatrixType::Lower) + { + const Complex *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'L'; + char dia = 'N'; + + Array<Complex> z (2 * nc); + Complex *pz = z.fortran_vec (); + Array<double> rz (nc); + double *prz = rz.fortran_vec (); + + F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, prz, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0.0; + } + else if (typ == MatrixType::Permuted_Lower) + (*current_liboctave_error_handler) + ("permuted triangular matrix not implemented"); + else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) + { + double anorm = -1.0; + ComplexMatrix atmp = *this; + Complex *tmp_data = atmp.fortran_vec (); + + if (typ == MatrixType::Hermitian) + { + octave_idx_type info = 0; + char job = 'L'; + anorm = atmp.abs().sum(). + row(static_cast<octave_idx_type>(0)).max(); + + F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, + tmp_data, nr, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + { + rcon = 0.0; + + mattype.mark_as_unsymmetric (); + typ = MatrixType::Full; + } + else + { + Array<Complex> z (2 * nc); + Complex *pz = z.fortran_vec (); + Array<double> rz (nc); + double *prz = rz.fortran_vec (); + + F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, tmp_data, nr, anorm, + rcon, pz, prz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0.0; + } + } + + + if (typ == MatrixType::Full) + { + octave_idx_type info = 0; + + Array<octave_idx_type> ipvt (nr); + octave_idx_type *pipvt = ipvt.fortran_vec (); + + if(anorm < 0.) + anorm = atmp.abs().sum(). + row(static_cast<octave_idx_type>(0)).max(); + + Array<Complex> z (2 * nc); + Complex *pz = z.fortran_vec (); + Array<double> rz (2 * nc); + double *prz = rz.fortran_vec (); + + F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info)); + + if (info != 0) + { + rcon = 0.0; + mattype.mark_as_rectangular (); + } + else + { + char job = '1'; + F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcon, pz, prz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0.0; + } + } + } + else + rcon = 0.0; + } + + return rcon; +} + ComplexMatrix ComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { @@ -1682,7 +1849,7 @@ typ == MatrixType::Upper) { octave_idx_type b_nc = b.cols (); - rcond = 1.; + rcon = 1.; info = 0; if (typ == MatrixType::Permuted_Upper) @@ -1708,7 +1875,7 @@ F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), F77_CONST_CHAR_ARG2 (&uplo, 1), F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcond, + nr, tmp_data, nr, rcon, pz, prz, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) @@ -1717,18 +1884,18 @@ if (info != 0) info = -2; - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile double rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) { info = -2; if (sing_handler) - sing_handler (rcond); + sing_handler (rcon); else (*current_liboctave_error_handler) ("matrix singular to machine precision, rcond = %g", - rcond); + rcon); } } @@ -1761,7 +1928,7 @@ ComplexMatrix ComplexMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { @@ -1783,7 +1950,7 @@ typ == MatrixType::Lower) { octave_idx_type b_nc = b.cols (); - rcond = 1.; + rcon = 1.; info = 0; if (typ == MatrixType::Permuted_Lower) @@ -1809,7 +1976,7 @@ F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), F77_CONST_CHAR_ARG2 (&uplo, 1), F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcond, + nr, tmp_data, nr, rcon, pz, prz, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) @@ -1818,18 +1985,18 @@ if (info != 0) info = -2; - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile double rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) { info = -2; if (sing_handler) - sing_handler (rcond); + sing_handler (rcon); else (*current_liboctave_error_handler) ("matrix singular to machine precision, rcond = %g", - rcond); + rcon); } } @@ -1862,7 +2029,7 @@ ComplexMatrix ComplexMatrix::fsolve (MatrixType &mattype, const ComplexMatrix& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { @@ -1897,7 +2064,7 @@ F77_CHAR_ARG_LEN (1))); // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + rcon = 0.0; if (info != 0) { info = -2; @@ -1916,24 +2083,24 @@ F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), nr, tmp_data, nr, anorm, - rcond, pz, prz, info + rcon, pz, prz, info F77_CHAR_ARG_LEN (1))); if (info != 0) info = -2; - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile double rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) { info = -2; if (sing_handler) - sing_handler (rcond); + sing_handler (rcon); else (*current_liboctave_error_handler) ("matrix singular to machine precision, rcond = %g", - rcond); + rcon); } } @@ -1979,13 +2146,13 @@ F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info)); // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + rcon = 0.0; if (info != 0) { info = -2; if (sing_handler) - sing_handler (rcond); + sing_handler (rcon); else (*current_liboctave_error_handler) ("matrix singular to machine precision"); @@ -2001,24 +2168,24 @@ char job = '1'; F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, tmp_data, nr, anorm, - rcond, pz, prz, info + rcon, pz, prz, info F77_CHAR_ARG_LEN (1))); if (info != 0) info = -2; - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile double rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) { info = -2; if (sing_handler) - sing_handler (rcond); + sing_handler (rcon); else (*current_liboctave_error_handler) ("matrix singular to machine precision, rcond = %g", - rcond); + rcon); } } @@ -2048,60 +2215,60 @@ ComplexMatrix::solve (MatrixType &typ, const Matrix& b) const { octave_idx_type info; - double rcond; - return solve (typ, b, info, rcond, 0); + double rcon; + return solve (typ, b, info, rcon, 0); } ComplexMatrix ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info) const { - double rcond; - return solve (typ, b, info, rcond, 0); + double rcon; + return solve (typ, b, info, rcon, 0); } ComplexMatrix ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, - double& rcond) const + double& rcon) const { - return solve (typ, b, info, rcond, 0); + return solve (typ, b, info, rcon, 0); } ComplexMatrix ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler, + double& rcon, solve_singularity_handler sing_handler, bool singular_fallback) const { ComplexMatrix tmp (b); - return solve (typ, tmp, info, rcond, sing_handler, singular_fallback); + return solve (typ, tmp, info, rcon, sing_handler, singular_fallback); } ComplexMatrix ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b) const { octave_idx_type info; - double rcond; - return solve (typ, b, info, rcond, 0); + double rcon; + return solve (typ, b, info, rcon, 0); } ComplexMatrix ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info) const { - double rcond; - return solve (typ, b, info, rcond, 0); + double rcon; + return solve (typ, b, info, rcon, 0); } ComplexMatrix ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b, - octave_idx_type& info, double& rcond) const + octave_idx_type& info, double& rcon) const { - return solve (typ, b, info, rcond, 0); + return solve (typ, b, info, rcon, 0); } ComplexMatrix ComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, bool singular_fallback) const { @@ -2113,11 +2280,11 @@ // Only calculate the condition number for LU/Cholesky if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) - retval = utsolve (mattype, b, info, rcond, sing_handler, false); + retval = utsolve (mattype, b, info, rcon, sing_handler, false); else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) - retval = ltsolve (mattype, b, info, rcond, sing_handler, false); + retval = ltsolve (mattype, b, info, rcon, sing_handler, false); else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) - retval = fsolve (mattype, b, info, rcond, sing_handler, true); + retval = fsolve (mattype, b, info, rcon, sing_handler, true); else if (typ != MatrixType::Rectangular) { (*current_liboctave_error_handler) ("unknown matrix type"); @@ -2128,7 +2295,7 @@ if (singular_fallback && mattype.type () == MatrixType::Rectangular) { octave_idx_type rank; - retval = lssolve (b, info, rank, rcond); + retval = lssolve (b, info, rank, rcon); } return retval; @@ -2138,183 +2305,183 @@ ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b) const { octave_idx_type info; - double rcond; - return solve (typ, ComplexColumnVector (b), info, rcond, 0); + double rcon; + return solve (typ, ComplexColumnVector (b), info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info) const { - double rcond; - return solve (typ, ComplexColumnVector (b), info, rcond, 0); + double rcon; + return solve (typ, ComplexColumnVector (b), info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b, - octave_idx_type& info, double& rcond) const + octave_idx_type& info, double& rcon) const { - return solve (typ, ComplexColumnVector (b), info, rcond, 0); + return solve (typ, ComplexColumnVector (b), info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const { - return solve (typ, ComplexColumnVector (b), info, rcond, sing_handler); + return solve (typ, ComplexColumnVector (b), info, rcon, sing_handler); } ComplexColumnVector ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b) const { octave_idx_type info; - double rcond; - return solve (typ, b, info, rcond, 0); + double rcon; + return solve (typ, b, info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b, octave_idx_type& info) const { - double rcond; - return solve (typ, b, info, rcond, 0); + double rcon; + return solve (typ, b, info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b, - octave_idx_type& info, double& rcond) const + octave_idx_type& info, double& rcon) const { - return solve (typ, b, info, rcond, 0); + return solve (typ, b, info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const { ComplexMatrix tmp (b); - return solve (typ, tmp, info, rcond, sing_handler).column(static_cast<octave_idx_type> (0)); + return solve (typ, tmp, info, rcon, sing_handler).column(static_cast<octave_idx_type> (0)); } ComplexMatrix ComplexMatrix::solve (const Matrix& b) const { octave_idx_type info; - double rcond; - return solve (b, info, rcond, 0); + double rcon; + return solve (b, info, rcon, 0); } ComplexMatrix ComplexMatrix::solve (const Matrix& b, octave_idx_type& info) const { - double rcond; - return solve (b, info, rcond, 0); + double rcon; + return solve (b, info, rcon, 0); } ComplexMatrix -ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcond) const +ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcon) const { - return solve (b, info, rcond, 0); + return solve (b, info, rcon, 0); } ComplexMatrix -ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcond, +ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const { ComplexMatrix tmp (b); - return solve (tmp, info, rcond, sing_handler); + return solve (tmp, info, rcon, sing_handler); } ComplexMatrix ComplexMatrix::solve (const ComplexMatrix& b) const { octave_idx_type info; - double rcond; - return solve (b, info, rcond, 0); + double rcon; + return solve (b, info, rcon, 0); } ComplexMatrix ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info) const { - double rcond; - return solve (b, info, rcond, 0); + double rcon; + return solve (b, info, rcon, 0); } ComplexMatrix -ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcond) const +ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon) const { - return solve (b, info, rcond, 0); + return solve (b, info, rcon, 0); } ComplexMatrix -ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcond, +ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcond, sing_handler); + return solve (mattype, b, info, rcon, sing_handler); } ComplexColumnVector ComplexMatrix::solve (const ColumnVector& b) const { octave_idx_type info; - double rcond; - return solve (ComplexColumnVector (b), info, rcond, 0); + double rcon; + return solve (ComplexColumnVector (b), info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info) const { - double rcond; - return solve (ComplexColumnVector (b), info, rcond, 0); + double rcon; + return solve (ComplexColumnVector (b), info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, - double& rcond) const + double& rcon) const { - return solve (ComplexColumnVector (b), info, rcond, 0); + return solve (ComplexColumnVector (b), info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, - double& rcond, + double& rcon, solve_singularity_handler sing_handler) const { - return solve (ComplexColumnVector (b), info, rcond, sing_handler); + return solve (ComplexColumnVector (b), info, rcon, sing_handler); } ComplexColumnVector ComplexMatrix::solve (const ComplexColumnVector& b) const { octave_idx_type info; - double rcond; - return solve (b, info, rcond, 0); + double rcon; + return solve (b, info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info) const { - double rcond; - return solve (b, info, rcond, 0); + double rcon; + return solve (b, info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, - double& rcond) const + double& rcon) const { - return solve (b, info, rcond, 0); + return solve (b, info, rcon, 0); } ComplexColumnVector ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, - double& rcond, + double& rcon, solve_singularity_handler sing_handler) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcond, sing_handler); + return solve (mattype, b, info, rcon, sing_handler); } ComplexMatrix @@ -2322,31 +2489,31 @@ { octave_idx_type info; octave_idx_type rank; - double rcond; - return lssolve (ComplexMatrix (b), info, rank, rcond); + double rcon; + return lssolve (ComplexMatrix (b), info, rank, rcon); } ComplexMatrix ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info) const { octave_idx_type rank; - double rcond; - return lssolve (ComplexMatrix (b), info, rank, rcond); + double rcon; + return lssolve (ComplexMatrix (b), info, rank, rcon); } ComplexMatrix ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type& rank) const { - double rcond; - return lssolve (ComplexMatrix (b), info, rank, rcond); + double rcon; + return lssolve (ComplexMatrix (b), info, rank, rcon); } ComplexMatrix ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info, - octave_idx_type& rank, double& rcond) const + octave_idx_type& rank, double& rcon) const { - return lssolve (ComplexMatrix (b), info, rank, rcond); + return lssolve (ComplexMatrix (b), info, rank, rcon); } ComplexMatrix @@ -2354,29 +2521,29 @@ { octave_idx_type info; octave_idx_type rank; - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } ComplexMatrix ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info) const { octave_idx_type rank; - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } ComplexMatrix ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, octave_idx_type& rank) const { - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } ComplexMatrix ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank, double& rcond) const + octave_idx_type& rank, double& rcon) const { ComplexMatrix retval; @@ -2394,7 +2561,7 @@ { volatile octave_idx_type minmn = (m < n ? m : n); octave_idx_type maxmn = m > n ? m : n; - rcond = -1.0; + rcon = -1.0; if (m != n) { @@ -2461,7 +2628,7 @@ octave_idx_type* piwork = iwork.fortran_vec (); F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, - ps, rcond, rank, work.fortran_vec (), + ps, rcon, rank, work.fortran_vec (), lwork, prwork, piwork, info)); // The workspace query is broken in at least LAPACK 3.0.0 @@ -2498,19 +2665,19 @@ work.resize (lwork); F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcond, rank, + maxmn, ps, rcon, rank, work.fortran_vec (), lwork, prwork, piwork, info)); if (rank < minmn) (*current_liboctave_warning_handler) ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e", - m, n, rank, rcond); + m, n, rank, rcon); if (s.elem (0) == 0.0) - rcond = 0.0; + rcon = 0.0; else - rcond = s.elem (minmn - 1) / s.elem (0); + rcon = s.elem (minmn - 1) / s.elem (0); retval.resize (n, nrhs); } @@ -2523,31 +2690,31 @@ { octave_idx_type info; octave_idx_type rank; - double rcond; - return lssolve (ComplexColumnVector (b), info, rank, rcond); + double rcon; + return lssolve (ComplexColumnVector (b), info, rank, rcon); } ComplexColumnVector ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info) const { octave_idx_type rank; - double rcond; - return lssolve (ComplexColumnVector (b), info, rank, rcond); + double rcon; + return lssolve (ComplexColumnVector (b), info, rank, rcon); } ComplexColumnVector ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const { - double rcond; - return lssolve (ComplexColumnVector (b), info, rank, rcond); + double rcon; + return lssolve (ComplexColumnVector (b), info, rank, rcon); } ComplexColumnVector ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, double& rcond) const + octave_idx_type& rank, double& rcon) const { - return lssolve (ComplexColumnVector (b), info, rank, rcond); + return lssolve (ComplexColumnVector (b), info, rank, rcon); } ComplexColumnVector @@ -2555,30 +2722,30 @@ { octave_idx_type info; octave_idx_type rank; - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } ComplexColumnVector ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info) const { octave_idx_type rank; - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } ComplexColumnVector ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const { - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } ComplexColumnVector ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, double& rcond) const + octave_idx_type& rank, double& rcon) const { ComplexColumnVector retval; @@ -2596,7 +2763,7 @@ { volatile octave_idx_type minmn = (m < n ? m : n); octave_idx_type maxmn = m > n ? m : n; - rcond = -1.0; + rcon = -1.0; if (m != n) { @@ -2655,7 +2822,7 @@ octave_idx_type* piwork = iwork.fortran_vec (); F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, - ps, rcond, rank, work.fortran_vec (), + ps, rcon, rank, work.fortran_vec (), lwork, prwork, piwork, info)); lwork = static_cast<octave_idx_type> (std::real (work(0))); @@ -2664,7 +2831,7 @@ iwork.resize (iwork(0)); F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcond, rank, + maxmn, ps, rcon, rank, work.fortran_vec (), lwork, prwork, piwork, info)); @@ -2673,12 +2840,12 @@ if (rank < minmn) (*current_liboctave_warning_handler) ("zgelsd: rank deficient %dx%d matrix, rank = %d, tol = %e", - m, n, rank, rcond); + m, n, rank, rcon); if (s.elem (0) == 0.0) - rcond = 0.0; + rcon = 0.0; else - rcond = s.elem (minmn - 1) / s.elem (0); + rcon = s.elem (minmn - 1) / s.elem (0); retval.resize (n, nrhs); } @@ -2702,11 +2869,11 @@ }; static void -solve_singularity_warning (double rcond) +solve_singularity_warning (double rcon) { (*current_liboctave_warning_handler) ("singular matrix encountered in expm calculation, rcond = %g", - rcond); + rcon); } ComplexMatrix @@ -2833,8 +3000,8 @@ // Compute pade approximation = inverse (dpp) * npp. - double rcond; - retval = dpp.solve (npp, info, rcond, solve_singularity_warning); + double rcon; + retval = dpp.solve (npp, info, rcon, solve_singularity_warning); if (info < 0) return retval;
--- a/liboctave/CMatrix.h Tue May 13 21:12:12 2008 +0200 +++ b/liboctave/CMatrix.h Wed May 14 18:09:56 2008 +0200 @@ -38,7 +38,7 @@ { public: - typedef void (*solve_singularity_handler) (double rcond); + typedef void (*solve_singularity_handler) (double rcon); ComplexMatrix (void) : MArray2<Complex> () { } @@ -142,21 +142,21 @@ private: ComplexMatrix tinverse (MatrixType &mattype, octave_idx_type& info, - double& rcond, int force, int calc_cond) const; + double& rcon, int force, int calc_cond) const; ComplexMatrix finverse (MatrixType &mattype, octave_idx_type& info, - double& rcond, int force, int calc_cond) const; + double& rcon, int force, int calc_cond) const; public: ComplexMatrix inverse (void) const; ComplexMatrix inverse (octave_idx_type& info) const; - ComplexMatrix inverse (octave_idx_type& info, double& rcond, int force = 0, + ComplexMatrix inverse (octave_idx_type& info, double& rcon, int force = 0, int calc_cond = 1) const; ComplexMatrix inverse (MatrixType &mattype) const; ComplexMatrix inverse (MatrixType &mattype, octave_idx_type& info) const; ComplexMatrix inverse (MatrixType &mattype, octave_idx_type& info, - double& rcond, int force = 0, + double& rcon, int force = 0, int calc_cond = 1) const; ComplexMatrix pseudo_inverse (double tol = 0.0) const; @@ -169,24 +169,27 @@ ComplexDET determinant (void) const; ComplexDET determinant (octave_idx_type& info) const; - ComplexDET determinant (octave_idx_type& info, double& rcond, int calc_cond = 1) const; + ComplexDET determinant (octave_idx_type& info, double& rcon, int calc_cond = 1) const; + + double rcond (void) const; + double rcond (MatrixType &mattype) const; private: // Upper triangular matrix solvers ComplexMatrix utsolve (MatrixType &typ, const ComplexMatrix& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; // Lower triangular matrix solvers ComplexMatrix ltsolve (MatrixType &typ, const ComplexMatrix& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; // Full matrix solvers (umfpack/cholesky) ComplexMatrix fsolve (MatrixType &typ, const ComplexMatrix& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; @@ -196,18 +199,18 @@ ComplexMatrix solve (MatrixType &typ, const Matrix& b, octave_idx_type& info) const; ComplexMatrix solve (MatrixType &typ, const Matrix& b, - octave_idx_type& info, double& rcond) const; + octave_idx_type& info, double& rcon) const; ComplexMatrix solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler, + double& rcon, solve_singularity_handler sing_handler, bool singular_fallback = true) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b, - octave_idx_type& info, double& rcond) const; + octave_idx_type& info, double& rcon) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, bool singular_fallback = true) const; @@ -215,9 +218,9 @@ ComplexColumnVector solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info) const; ComplexColumnVector solve (MatrixType &typ, const ColumnVector& b, - octave_idx_type& info, double& rcond) const; + octave_idx_type& info, double& rcon) const; ComplexColumnVector solve (MatrixType &typ, const ColumnVector& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const; ComplexColumnVector solve (MatrixType &typ, @@ -225,37 +228,37 @@ ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b, octave_idx_type& info) const; ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b, - octave_idx_type& info, double& rcond) const; + octave_idx_type& info, double& rcon) const; ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const; // 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, + ComplexMatrix solve (const Matrix& b, octave_idx_type& info, double& rcon) const; + ComplexMatrix solve (const Matrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const; ComplexMatrix solve (const ComplexMatrix& b) const; 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, + ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon) const; + ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon, 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, + double& rcon) const; + ComplexColumnVector solve (const ColumnVector& b, octave_idx_type& info, double& rcon, 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, - double& rcond) const; + double& rcon) const; ComplexColumnVector solve (const ComplexColumnVector& b, octave_idx_type& info, - double& rcond, + double& rcon, solve_singularity_handler sing_handler) const; ComplexMatrix lssolve (const Matrix& b) const; @@ -263,14 +266,14 @@ 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; + octave_idx_type& rank, double& rcon) 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; + octave_idx_type& rank, double& rcon) const; ComplexColumnVector lssolve (const ColumnVector& b) const; ComplexColumnVector lssolve (const ColumnVector& b, @@ -278,7 +281,7 @@ 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; + octave_idx_type& rank, double& rcon) const; ComplexColumnVector lssolve (const ComplexColumnVector& b) const; ComplexColumnVector lssolve (const ComplexColumnVector& b, @@ -288,7 +291,7 @@ octave_idx_type& rank) const; ComplexColumnVector lssolve (const ComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, double& rcond) const; + octave_idx_type& rank, double& rcon) const; ComplexMatrix expm (void) const;
--- a/liboctave/ChangeLog Tue May 13 21:12:12 2008 +0200 +++ b/liboctave/ChangeLog Wed May 14 18:09:56 2008 +0200 @@ -1,5 +1,26 @@ 2008-05-20 David Bateman <dbateman@free.fr> + * CMatrix.cc (double rcond): Replace with double rcon everywhere + to avoid shadowed variable warning + (double ComplexMatrix::rcond (void) const): New method for + reciprocal condition number calculation. + (double ComplexMatrix::rcond (MatrixType &mattype) const): ditto. + * CMatrix.h (double rcond): Replace with double rcon everywhere + to avoid shadowed variable warning + (double ComplexMatrix::rcond (void) const): New method for + reciprocal condition number calculation. + (double ComplexMatrix::rcond (MatrixType &mattype) const): ditto. + * dMatrix.cc (double rcond): Replace with double rcon everywhere + to avoid shadowed variable warning + (double Matrix::rcond (void) const): New method for + reciprocal condition number calculation. + (double Matrix::rcond (MatrixType &mattype) const): ditto. + * dMatrix.h (double rcond): Replace with double rcon everywhere + to avoid shadowed variable warning + (double Matrix::rcond (void) const): New method for + reciprocal condition number calculation. + (double Matrix::rcond (MatrixType &mattype) const): ditto. + * regex-match.cc, regex-match.h: New class for simple regular expression matching * Makefile.in (INCLUDES): Add regex-match.h here, and
--- a/liboctave/dMatrix.cc Tue May 13 21:12:12 2008 +0200 +++ b/liboctave/dMatrix.cc Wed May 14 18:09:56 2008 +0200 @@ -651,44 +651,44 @@ Matrix::inverse (void) const { octave_idx_type info; - double rcond; + double rcon; MatrixType mattype (*this); - return inverse (mattype, info, rcond, 0, 0); + return inverse (mattype, info, rcon, 0, 0); } Matrix Matrix::inverse (octave_idx_type& info) const { - double rcond; + double rcon; MatrixType mattype (*this); - return inverse (mattype, info, rcond, 0, 0); + return inverse (mattype, info, rcon, 0, 0); } Matrix -Matrix::inverse (octave_idx_type& info, double& rcond, int force, +Matrix::inverse (octave_idx_type& info, double& rcon, int force, int calc_cond) const { MatrixType mattype (*this); - return inverse (mattype, info, rcond, force, calc_cond); + return inverse (mattype, info, rcon, force, calc_cond); } Matrix Matrix::inverse (MatrixType& mattype) const { octave_idx_type info; - double rcond; - return inverse (mattype, info, rcond, 0, 0); + double rcon; + return inverse (mattype, info, rcon, 0, 0); } Matrix Matrix::inverse (MatrixType &mattype, octave_idx_type& info) const { - double rcond; - return inverse (mattype, info, rcond, 0, 0); + double rcon; + return inverse (mattype, info, rcon, 0, 0); } Matrix -Matrix::tinverse (MatrixType &mattype, octave_idx_type& info, double& rcond, +Matrix::tinverse (MatrixType &mattype, octave_idx_type& info, double& rcon, int force, int calc_cond) const { Matrix retval; @@ -713,7 +713,7 @@ F77_CHAR_ARG_LEN (1))); // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + rcon = 0.0; if (info != 0) info = -1; else if (calc_cond) @@ -727,7 +727,7 @@ F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&job, 1), F77_CONST_CHAR_ARG2 (&uplo, 1), F77_CONST_CHAR_ARG2 (&udiag, 1), - nr, tmp_data, nr, rcond, + nr, tmp_data, nr, rcon, work, iwork, dtrcon_info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) @@ -746,7 +746,7 @@ Matrix -Matrix::finverse (MatrixType &mattype, octave_idx_type& info, double& rcond, +Matrix::finverse (MatrixType &mattype, octave_idx_type& info, double& rcon, int force, int calc_cond) const { Matrix retval; @@ -786,7 +786,7 @@ F77_XFCN (dgetrf, DGETRF, (nc, nc, tmp_data, nr, pipvt, info)); // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + rcon = 0.0; if (info != 0) info = -1; else if (calc_cond) @@ -799,7 +799,7 @@ octave_idx_type *piz = iz.fortran_vec (); F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, tmp_data, nr, anorm, - rcond, pz, piz, dgecon_info + rcon, pz, piz, dgecon_info F77_CHAR_ARG_LEN (1))); if (dgecon_info != 0) @@ -827,7 +827,7 @@ } Matrix -Matrix::inverse (MatrixType &mattype, octave_idx_type& info, double& rcond, +Matrix::inverse (MatrixType &mattype, octave_idx_type& info, double& rcon, int force, int calc_cond) const { int typ = mattype.type (false); @@ -837,7 +837,7 @@ typ = mattype.type (*this); if (typ == MatrixType::Upper || typ == MatrixType::Lower) - ret = tinverse (mattype, info, rcond, force, calc_cond); + ret = tinverse (mattype, info, rcon, force, calc_cond); else { if (mattype.is_hermitian ()) @@ -846,9 +846,9 @@ if (info == 0) { if (calc_cond) - rcond = chol.rcond (); + rcon = chol.rcond (); else - rcond = 1.0; + rcon = 1.0; ret = chol.inverse (); } else @@ -856,9 +856,9 @@ } if (!mattype.is_hermitian ()) - ret = finverse(mattype, info, rcond, force, calc_cond); - - if ((mattype.is_hermitian () || calc_cond) && rcond == 0.) + ret = finverse(mattype, info, rcon, force, calc_cond); + + if ((mattype.is_hermitian () || calc_cond) && rcon == 0.) ret = Matrix (rows (), columns (), octave_Inf); } @@ -1216,19 +1216,19 @@ Matrix::determinant (void) const { octave_idx_type info; - double rcond; - return determinant (info, rcond, 0); + double rcon; + return determinant (info, rcon, 0); } DET Matrix::determinant (octave_idx_type& info) const { - double rcond; - return determinant (info, rcond, 0); + double rcon; + return determinant (info, rcon, 0); } DET -Matrix::determinant (octave_idx_type& info, double& rcond, int calc_cond) const +Matrix::determinant (octave_idx_type& info, double& rcon, int calc_cond) const { DET retval; @@ -1257,7 +1257,7 @@ F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + rcon = 0.0; if (info != 0) { info = -1; @@ -1276,7 +1276,7 @@ F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, tmp_data, nr, anorm, - rcond, pz, piz, info + rcon, pz, piz, info F77_CHAR_ARG_LEN (1))); } @@ -1321,9 +1321,174 @@ return retval; } +double +Matrix::rcond (void) const +{ + MatrixType mattype (*this); + return rcond (mattype); +} + +double +Matrix::rcond (MatrixType &mattype) const +{ + double rcon; + octave_idx_type nr = rows (); + octave_idx_type nc = cols (); + + if (nr != nc) + (*current_liboctave_error_handler) ("matrix must be square"); + else if (nr == 0 || nc == 0) + rcon = octave_Inf; + else + { + int typ = mattype.type (); + + if (typ == MatrixType::Unknown) + typ = mattype.type (*this); + + // Only calculate the condition number for LU/Cholesky + if (typ == MatrixType::Upper) + { + const double *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'U'; + char dia = 'N'; + + Array<double> z (3 * nc); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, piz, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0.0; + } + else if (typ == MatrixType::Permuted_Upper) + (*current_liboctave_error_handler) + ("permuted triangular matrix not implemented"); + else if (typ == MatrixType::Lower) + { + const double *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'L'; + char dia = 'N'; + + Array<double> z (3 * nc); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), + F77_CONST_CHAR_ARG2 (&uplo, 1), + F77_CONST_CHAR_ARG2 (&dia, 1), + nr, tmp_data, nr, rcon, + pz, piz, info + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1) + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0.0; + } + else if (typ == MatrixType::Permuted_Lower) + (*current_liboctave_error_handler) + ("permuted triangular matrix not implemented"); + else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) + { + double anorm = -1.0; + Matrix atmp = *this; + double *tmp_data = atmp.fortran_vec (); + + if (typ == MatrixType::Hermitian) + { + octave_idx_type info = 0; + char job = 'L'; + anorm = atmp.abs().sum(). + row(static_cast<octave_idx_type>(0)).max(); + + F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, + tmp_data, nr, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + { + rcon = 0.0; + mattype.mark_as_unsymmetric (); + typ = MatrixType::Full; + } + else + { + Array<double> z (3 * nc); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), + nr, tmp_data, nr, anorm, + rcon, pz, piz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0.0; + } + } + + if (typ == MatrixType::Full) + { + octave_idx_type info = 0; + + Array<octave_idx_type> ipvt (nr); + octave_idx_type *pipvt = ipvt.fortran_vec (); + + if(anorm < 0.) + anorm = atmp.abs().sum(). + row(static_cast<octave_idx_type>(0)).max(); + + Array<double> z (4 * nc); + double *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); + + if (info != 0) + { + rcon = 0.0; + mattype.mark_as_rectangular (); + } + else + { + char job = '1'; + F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), + nc, tmp_data, nr, anorm, + rcon, pz, piz, info + F77_CHAR_ARG_LEN (1))); + + if (info != 0) + rcon = 0.0; + } + } + } + else + rcon = 0.0; + } + + return rcon; +} + Matrix Matrix::utsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler, + double& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { Matrix retval; @@ -1344,7 +1509,7 @@ typ == MatrixType::Upper) { octave_idx_type b_nc = b.cols (); - rcond = 1.; + rcon = 1.; info = 0; if (typ == MatrixType::Permuted_Upper) @@ -1370,7 +1535,7 @@ F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), F77_CONST_CHAR_ARG2 (&uplo, 1), F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcond, + nr, tmp_data, nr, rcon, pz, piz, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) @@ -1379,18 +1544,18 @@ if (info != 0) info = -2; - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile double rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) { info = -2; if (sing_handler) - sing_handler (rcond); + sing_handler (rcon); else (*current_liboctave_error_handler) ("matrix singular to machine precision, rcond = %g", - rcond); + rcon); } } @@ -1423,7 +1588,7 @@ Matrix Matrix::ltsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler, + double& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { Matrix retval; @@ -1444,7 +1609,7 @@ typ == MatrixType::Lower) { octave_idx_type b_nc = b.cols (); - rcond = 1.; + rcon = 1.; info = 0; if (typ == MatrixType::Permuted_Lower) @@ -1470,7 +1635,7 @@ F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1), F77_CONST_CHAR_ARG2 (&uplo, 1), F77_CONST_CHAR_ARG2 (&dia, 1), - nr, tmp_data, nr, rcond, + nr, tmp_data, nr, rcon, pz, piz, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1) @@ -1479,18 +1644,18 @@ if (info != 0) info = -2; - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile double rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) { info = -2; if (sing_handler) - sing_handler (rcond); + sing_handler (rcon); else (*current_liboctave_error_handler) ("matrix singular to machine precision, rcond = %g", - rcond); + rcon); } } @@ -1523,7 +1688,7 @@ Matrix Matrix::fsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler, + double& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { Matrix retval; @@ -1556,7 +1721,7 @@ F77_CHAR_ARG_LEN (1))); // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + rcon = 0.0; if (info != 0) { info = -2; @@ -1575,24 +1740,24 @@ F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1), nr, tmp_data, nr, anorm, - rcond, pz, piz, info + rcon, pz, piz, info F77_CHAR_ARG_LEN (1))); if (info != 0) info = -2; - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile double rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) { info = -2; if (sing_handler) - sing_handler (rcond); + sing_handler (rcon); else (*current_liboctave_error_handler) ("matrix singular to machine precision, rcond = %g", - rcond); + rcon); } } @@ -1636,13 +1801,13 @@ F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info)); // Throw-away extra info LAPACK gives so as to not change output. - rcond = 0.0; + rcon = 0.0; if (info != 0) { info = -2; if (sing_handler) - sing_handler (rcond); + sing_handler (rcon); else (*current_liboctave_error_handler) ("matrix singular to machine precision"); @@ -1658,24 +1823,24 @@ char job = '1'; F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, tmp_data, nr, anorm, - rcond, pz, piz, info + rcon, pz, piz, info F77_CHAR_ARG_LEN (1))); if (info != 0) info = -2; - volatile double rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile double rcond_plus_one = rcon + 1.0; + + if (rcond_plus_one == 1.0 || xisnan (rcon)) { info = -2; if (sing_handler) - sing_handler (rcond); + sing_handler (rcon); else (*current_liboctave_error_handler) ("matrix singular to machine precision, rcond = %g", - rcond); + rcon); } } @@ -1707,20 +1872,20 @@ Matrix::solve (MatrixType &typ, const Matrix& b) const { octave_idx_type info; - double rcond; - return solve (typ, b, info, rcond, 0); + double rcon; + return solve (typ, b, info, rcon, 0); } Matrix Matrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, - double& rcond) const + double& rcon) const { - return solve (typ, b, info, rcond, 0); + return solve (typ, b, info, rcon, 0); } Matrix Matrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler, + double& rcon, solve_singularity_handler sing_handler, bool singular_fallback) const { Matrix retval; @@ -1731,11 +1896,11 @@ // Only calculate the condition number for LU/Cholesky if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) - retval = utsolve (mattype, b, info, rcond, sing_handler, false); + retval = utsolve (mattype, b, info, rcon, sing_handler, false); else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) - retval = ltsolve (mattype, b, info, rcond, sing_handler, false); + retval = ltsolve (mattype, b, info, rcon, sing_handler, false); else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) - retval = fsolve (mattype, b, info, rcond, sing_handler, true); + retval = fsolve (mattype, b, info, rcon, sing_handler, true); else if (typ != MatrixType::Rectangular) { (*current_liboctave_error_handler) ("unknown matrix type"); @@ -1746,7 +1911,7 @@ if (singular_fallback && mattype.type () == MatrixType::Rectangular) { octave_idx_type rank; - retval = lssolve (b, info, rank, rcond); + retval = lssolve (b, info, rank, rcon); } return retval; @@ -1769,49 +1934,49 @@ ComplexMatrix Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, - double& rcond) const + double& rcon) const { ComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcond); + return tmp.solve (typ, b, info, rcon); } ComplexMatrix Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler, + double& rcon, solve_singularity_handler sing_handler, bool singular_fallback) const { ComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcond, sing_handler, singular_fallback); + return tmp.solve (typ, b, info, rcon, sing_handler, singular_fallback); } ColumnVector Matrix::solve (MatrixType &typ, const ColumnVector& b) const { - octave_idx_type info; double rcond; - return solve (typ, b, info, rcond); + octave_idx_type info; double rcon; + return solve (typ, b, info, rcon); } ColumnVector Matrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info) const { - double rcond; - return solve (typ, b, info, rcond); + double rcon; + return solve (typ, b, info, rcon); } ColumnVector Matrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info, - double& rcond) const + double& rcon) const { - return solve (typ, b, info, rcond, 0); + return solve (typ, b, info, rcon, 0); } ColumnVector Matrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler) const + double& rcon, solve_singularity_handler sing_handler) const { Matrix tmp (b); - return solve (typ, tmp, info, rcond, sing_handler).column(static_cast<octave_idx_type> (0)); + return solve (typ, tmp, info, rcon, sing_handler).column(static_cast<octave_idx_type> (0)); } ComplexColumnVector @@ -1831,48 +1996,48 @@ ComplexColumnVector Matrix::solve (MatrixType &typ, const ComplexColumnVector& b, - octave_idx_type& info, double& rcond) const + octave_idx_type& info, double& rcon) const { ComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcond); + return tmp.solve (typ, b, info, rcon); } ComplexColumnVector Matrix::solve (MatrixType &typ, const ComplexColumnVector& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const { ComplexMatrix tmp (*this); - return tmp.solve(typ, b, info, rcond, sing_handler); + return tmp.solve(typ, b, info, rcon, sing_handler); } Matrix Matrix::solve (const Matrix& b) const { octave_idx_type info; - double rcond; - return solve (b, info, rcond, 0); + double rcon; + return solve (b, info, rcon, 0); } Matrix Matrix::solve (const Matrix& b, octave_idx_type& info) const { - double rcond; - return solve (b, info, rcond, 0); + double rcon; + return solve (b, info, rcon, 0); } Matrix -Matrix::solve (const Matrix& b, octave_idx_type& info, double& rcond) const +Matrix::solve (const Matrix& b, octave_idx_type& info, double& rcon) const { - return solve (b, info, rcond, 0); + return solve (b, info, rcon, 0); } Matrix Matrix::solve (const Matrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler) const + double& rcon, solve_singularity_handler sing_handler) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcond, sing_handler); + return solve (mattype, b, info, rcon, sing_handler); } ComplexMatrix @@ -1890,46 +2055,46 @@ } ComplexMatrix -Matrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcond) const +Matrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon) const { ComplexMatrix tmp (*this); - return tmp.solve (b, info, rcond); + return tmp.solve (b, info, rcon); } ComplexMatrix -Matrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcond, +Matrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const { ComplexMatrix tmp (*this); - return tmp.solve (b, info, rcond, sing_handler); + return tmp.solve (b, info, rcon, sing_handler); } ColumnVector Matrix::solve (const ColumnVector& b) const { - octave_idx_type info; double rcond; - return solve (b, info, rcond); + octave_idx_type info; double rcon; + return solve (b, info, rcon); } ColumnVector Matrix::solve (const ColumnVector& b, octave_idx_type& info) const { - double rcond; - return solve (b, info, rcond); + double rcon; + return solve (b, info, rcon); } ColumnVector -Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond) const +Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcon) const { - return solve (b, info, rcond, 0); + return solve (b, info, rcon, 0); } ColumnVector -Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond, +Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const { MatrixType mattype (*this); - return solve (mattype, b, info, rcond, sing_handler); + return solve (mattype, b, info, rcon, sing_handler); } ComplexColumnVector @@ -1947,18 +2112,18 @@ } ComplexColumnVector -Matrix::solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcond) const +Matrix::solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcon) const { ComplexMatrix tmp (*this); - return tmp.solve (b, info, rcond); + return tmp.solve (b, info, rcon); } ComplexColumnVector -Matrix::solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcond, +Matrix::solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const { ComplexMatrix tmp (*this); - return tmp.solve (b, info, rcond, sing_handler); + return tmp.solve (b, info, rcon, sing_handler); } Matrix @@ -1966,29 +2131,29 @@ { octave_idx_type info; octave_idx_type rank; - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } Matrix Matrix::lssolve (const Matrix& b, octave_idx_type& info) const { octave_idx_type rank; - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } Matrix Matrix::lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type& rank) const { - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } Matrix Matrix::lssolve (const Matrix& b, octave_idx_type& info, - octave_idx_type& rank, double &rcond) const + octave_idx_type& rank, double &rcon) const { Matrix retval; @@ -2006,7 +2171,7 @@ { volatile octave_idx_type minmn = (m < n ? m : n); octave_idx_type maxmn = m > n ? m : n; - rcond = -1.0; + rcon = -1.0; if (m != n) { retval = Matrix (maxmn, nrhs, 0.0); @@ -2064,7 +2229,7 @@ octave_idx_type* piwork = iwork.fortran_vec (); F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, - ps, rcond, rank, work.fortran_vec (), + ps, rcon, rank, work.fortran_vec (), lwork, piwork, info)); // The workspace query is broken in at least LAPACK 3.0.0 @@ -2108,7 +2273,7 @@ work.resize (lwork); F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcond, rank, + maxmn, ps, rcon, rank, work.fortran_vec (), lwork, piwork, info)); @@ -2116,9 +2281,9 @@ (*current_liboctave_warning_handler) ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank); if (s.elem (0) == 0.0) - rcond = 0.0; + rcon = 0.0; else - rcond = s.elem (minmn - 1) / s.elem (0); + rcon = s.elem (minmn - 1) / s.elem (0); retval.resize (n, nrhs); } @@ -2132,8 +2297,8 @@ ComplexMatrix tmp (*this); octave_idx_type info; octave_idx_type rank; - double rcond; - return tmp.lssolve (b, info, rank, rcond); + double rcon; + return tmp.lssolve (b, info, rank, rcon); } ComplexMatrix @@ -2141,8 +2306,8 @@ { ComplexMatrix tmp (*this); octave_idx_type rank; - double rcond; - return tmp.lssolve (b, info, rank, rcond); + double rcon; + return tmp.lssolve (b, info, rank, rcon); } ComplexMatrix @@ -2150,16 +2315,16 @@ octave_idx_type& rank) const { ComplexMatrix tmp (*this); - double rcond; - return tmp.lssolve (b, info, rank, rcond); + double rcon; + return tmp.lssolve (b, info, rank, rcon); } ComplexMatrix Matrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank, double& rcond) const + octave_idx_type& rank, double& rcon) const { ComplexMatrix tmp (*this); - return tmp.lssolve (b, info, rank, rcond); + return tmp.lssolve (b, info, rank, rcon); } ColumnVector @@ -2167,29 +2332,29 @@ { octave_idx_type info; octave_idx_type rank; - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } ColumnVector Matrix::lssolve (const ColumnVector& b, octave_idx_type& info) const { octave_idx_type rank; - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } ColumnVector Matrix::lssolve (const ColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const { - double rcond; - return lssolve (b, info, rank, rcond); + double rcon; + return lssolve (b, info, rank, rcon); } ColumnVector Matrix::lssolve (const ColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, double &rcond) const + octave_idx_type& rank, double &rcon) const { ColumnVector retval; @@ -2207,7 +2372,7 @@ { volatile octave_idx_type minmn = (m < n ? m : n); octave_idx_type maxmn = m > n ? m : n; - rcond = -1.0; + rcon = -1.0; if (m != n) { @@ -2258,14 +2423,14 @@ octave_idx_type* piwork = iwork.fortran_vec (); F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, - ps, rcond, rank, work.fortran_vec (), + ps, rcon, rank, work.fortran_vec (), lwork, piwork, info)); lwork = static_cast<octave_idx_type> (work(0)); work.resize (lwork); F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcond, rank, + maxmn, ps, rcon, rank, work.fortran_vec (), lwork, piwork, info)); @@ -2275,9 +2440,9 @@ (*current_liboctave_warning_handler) ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank); if (s.elem (0) == 0.0) - rcond = 0.0; + rcon = 0.0; else - rcond = s.elem (minmn - 1) / s.elem (0); + rcon = s.elem (minmn - 1) / s.elem (0); } retval.resize (n, nrhs); @@ -2292,8 +2457,8 @@ ComplexMatrix tmp (*this); octave_idx_type info; octave_idx_type rank; - double rcond; - return tmp.lssolve (b, info, rank, rcond); + double rcon; + return tmp.lssolve (b, info, rank, rcon); } ComplexColumnVector @@ -2301,8 +2466,8 @@ { ComplexMatrix tmp (*this); octave_idx_type rank; - double rcond; - return tmp.lssolve (b, info, rank, rcond); + double rcon; + return tmp.lssolve (b, info, rank, rcon); } ComplexColumnVector @@ -2310,16 +2475,16 @@ octave_idx_type& rank) const { ComplexMatrix tmp (*this); - double rcond; - return tmp.lssolve (b, info, rank, rcond); + double rcon; + return tmp.lssolve (b, info, rank, rcon); } ComplexColumnVector Matrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, double &rcond) const + octave_idx_type& rank, double &rcon) const { ComplexMatrix tmp (*this); - return tmp.lssolve (b, info, rank, rcond); + return tmp.lssolve (b, info, rank, rcon); } // Constants for matrix exponential calculation. @@ -2337,11 +2502,11 @@ }; static void -solve_singularity_warning (double rcond) +solve_singularity_warning (double rcon) { (*current_liboctave_warning_handler) ("singular matrix encountered in expm calculation, rcond = %g", - rcond); + rcon); } Matrix @@ -2466,8 +2631,8 @@ // Compute pade approximation = inverse (dpp) * npp. - double rcond; - retval = dpp.solve (npp, info, rcond, solve_singularity_warning); + double rcon; + retval = dpp.solve (npp, info, rcon, solve_singularity_warning); if (info < 0) return retval;
--- a/liboctave/dMatrix.h Tue May 13 21:12:12 2008 +0200 +++ b/liboctave/dMatrix.h Wed May 14 18:09:56 2008 +0200 @@ -37,7 +37,7 @@ { public: - typedef void (*solve_singularity_handler) (double rcond); + typedef void (*solve_singularity_handler) (double rcon); Matrix (void) : MArray2<double> () { } @@ -112,21 +112,21 @@ ColumnVector column (octave_idx_type i) const; private: - Matrix tinverse (MatrixType &mattype, octave_idx_type& info, double& rcond, + Matrix tinverse (MatrixType &mattype, octave_idx_type& info, double& rcon, int force, int calc_cond) const; - Matrix finverse (MatrixType &mattype, octave_idx_type& info, double& rcond, + Matrix finverse (MatrixType &mattype, octave_idx_type& info, double& rcon, int force, int calc_cond) const; public: Matrix inverse (void) const; Matrix inverse (octave_idx_type& info) const; - Matrix inverse (octave_idx_type& info, double& rcond, int force = 0, + Matrix inverse (octave_idx_type& info, double& rcon, int force = 0, int calc_cond = 1) const; Matrix inverse (MatrixType &mattype) const; Matrix inverse (MatrixType &mattype, octave_idx_type& info) const; - Matrix inverse (MatrixType &mattype, octave_idx_type& info, double& rcond, + Matrix inverse (MatrixType &mattype, octave_idx_type& info, double& rcon, int force = 0, int calc_cond = 1) const; Matrix pseudo_inverse (double tol = 0.0) const; @@ -139,22 +139,25 @@ DET determinant (void) const; DET determinant (octave_idx_type& info) const; - DET determinant (octave_idx_type& info, double& rcond, int calc_cond = 1) const; + DET determinant (octave_idx_type& info, double& rcon, int calc_cond = 1) const; + + double rcond (void) const; + double rcond (MatrixType &mattype) const; private: // Upper triangular matrix solvers Matrix utsolve (MatrixType &typ, const Matrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler, + double& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; // Lower triangular matrix solvers Matrix ltsolve (MatrixType &typ, const Matrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler, + double& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; // Full matrix solvers (lu/cholesky) Matrix fsolve (MatrixType &typ, const Matrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler, + double& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; public: @@ -162,18 +165,18 @@ Matrix solve (MatrixType &typ, const Matrix& b) const; Matrix solve (MatrixType &typ, const Matrix& b, octave_idx_type& info) const; Matrix solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, - double& rcond) const; + double& rcon) const; Matrix solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, - double& rcond, solve_singularity_handler sing_handler, + double& rcon, solve_singularity_handler sing_handler, bool singular_fallback = true) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b, - octave_idx_type& info, double& rcond) const; + octave_idx_type& info, double& rcon) const; ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler, bool singular_fallback = true) const; @@ -181,9 +184,9 @@ ColumnVector solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info) const; ColumnVector solve (MatrixType &typ, const ColumnVector& b, - octave_idx_type& info, double& rcond) const; + octave_idx_type& info, double& rcon) const; ColumnVector solve (MatrixType &typ, const ColumnVector& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const; ComplexColumnVector solve (MatrixType &typ, @@ -191,36 +194,36 @@ ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b, octave_idx_type& info) const; ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b, - octave_idx_type& info, double& rcond) const; + octave_idx_type& info, double& rcon) const; ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b, - octave_idx_type& info, double& rcond, + octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const; // Generic interface to solver with probing of type Matrix solve (const Matrix& b) const; Matrix solve (const Matrix& b, octave_idx_type& info) const; - Matrix solve (const Matrix& b, octave_idx_type& info, double& rcond) const; - Matrix solve (const Matrix& b, octave_idx_type& info, double& rcond, + Matrix solve (const Matrix& b, octave_idx_type& info, double& rcon) const; + Matrix solve (const Matrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const; ComplexMatrix solve (const ComplexMatrix& b) const; 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, + ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon) const; + ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon, solve_singularity_handler sing_handler) const; ColumnVector solve (const ColumnVector& b) const; ColumnVector solve (const ColumnVector& b, octave_idx_type& info) const; - ColumnVector solve (const ColumnVector& b, octave_idx_type& info, double& rcond) const; - ColumnVector solve (const ColumnVector& b, octave_idx_type& info, double& rcond, + ColumnVector solve (const ColumnVector& b, octave_idx_type& info, double& rcon) const; + ColumnVector solve (const ColumnVector& b, octave_idx_type& info, double& rcon, 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, - double& rcond) const; + double& rcon) const; ComplexColumnVector solve (const ComplexColumnVector& b, octave_idx_type& info, - double& rcond, + double& rcon, solve_singularity_handler sing_handler) const; // Singular solvers @@ -229,21 +232,21 @@ 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; + octave_idx_type& rank, double& rcon) 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; + octave_idx_type& rank, double &rcon) 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, double& rcond) const; + octave_idx_type& rank, double& rcon) const; ComplexColumnVector lssolve (const ComplexColumnVector& b) const; ComplexColumnVector lssolve (const ComplexColumnVector& b, @@ -253,7 +256,7 @@ octave_idx_type& rank) const; ComplexColumnVector lssolve (const ComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, double& rcond) const; + octave_idx_type& rank, double& rcon) const; Matrix expm (void) const;
--- a/src/ChangeLog Tue May 13 21:12:12 2008 +0200 +++ b/src/ChangeLog Wed May 14 18:09:56 2008 +0200 @@ -1,5 +1,8 @@ 2008-05-20 David Bateman <dbateman@free.fr> + * DLD-FUNCTIONS/rcond.cc: New function. + * Makefile.in (DLD_XSRC): Add it here. + * debug.cc (Fdbstop): If no line specified assume line 1. (Fdbstep, Fdbcont, Fdbnext): Move debugging functions to normal commands.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/rcond.cc Wed May 14 18:09:56 2008 +0200 @@ -0,0 +1,70 @@ +/* + +Copyright (C) 2008 David Bateman + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (rcond, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{c} =} rcond (@var{a})\n\ +Compute the 1-norm estimate of the reciprocal condition as returned\n\ +by LAPACK. If the matrix is well-conditioned then @var{c} will be near\n\ +1 and if the matrix is poorly conditioned it will be close to zero.\n\ +\n\ +The matrix @var{a} must not be sparse. If the matrix is sparse then\n\ +@code{condest (@var{a})} or @code{rcond (full (@var{a}))} should be used\n\ +instead.\n\ +@seealso{inv, mldivide}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin != 1) + print_usage (); + else if (args(0).is_sparse_type ()) + error ("rcond: for sparse matrices use 'rcond (full (a))' or 'condest (a)' instead"); + else if (args(0).is_complex_type ()) + { + ComplexMatrix m = args(0).complex_matrix_value (); + MatrixType mattyp; + retval = m.rcond (mattyp); + args(0).matrix_type (mattyp); + } + else + { + Matrix m = args(0).matrix_value (); + MatrixType mattyp; + retval = m.rcond (mattyp); + args(0).matrix_type (mattyp); + } + + return retval; +}
--- a/src/Makefile.in Tue May 13 21:12:12 2008 +0200 +++ b/src/Makefile.in Wed May 14 18:09:56 2008 +0200 @@ -69,7 +69,7 @@ gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \ givens.cc hess.cc hex2num.cc inv.cc kron.cc lookup.cc lsode.cc \ lu.cc luinc.cc matrix_type.cc max.cc md5sum.cc pinv.cc qr.cc \ - quad.cc qz.cc rand.cc regexp.cc schur.cc sparse.cc \ + quad.cc qz.cc rand.cc rcond.cc regexp.cc schur.cc sparse.cc \ spparms.cc sqrtm.cc svd.cc syl.cc symrcm.cc symbfact.cc \ time.cc tsearch.cc typecast.cc \ urlwrite.cc __contourc__.cc __delaunayn__.cc __dsearchn__.cc \