Mercurial > octave-nkf
changeset 7797:f42c6f8d6d8e
Extend rcond function to single precision types
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Wed, 14 May 2008 23:28:41 +0200 |
| parents | 762801c50b21 |
| children | 42c42c640108 |
| files | liboctave/ChangeLog liboctave/fCMatrix.cc liboctave/fCMatrix.h liboctave/fMatrix.cc liboctave/fMatrix.h src/ChangeLog src/DLD-FUNCTIONS/rcond.cc |
| diffstat | 7 files changed, 766 insertions(+), 388 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog Tue May 13 12:40:23 2008 +0200 +++ b/liboctave/ChangeLog Wed May 14 23:28:41 2008 +0200 @@ -1,5 +1,26 @@ 2008-05-21 David Bateman <dbateman@free.fr> + * fCMatrix.cc (float rcond): Replace with float rcon everywhere + to avoid shadowed variable warning + (float ComplexMatrix::rcond (void) const): New method for + reciprocal condition number calculation. + (float ComplexMatrix::rcond (MatrixType &mattype) const): ditto. + * fCMatrix.h (float rcond): Replace with float rcon everywhere + to avoid shadowed variable warning + (float ComplexMatrix::rcond (void) const): New method for + reciprocal condition number calculation. + (float ComplexMatrix::rcond (MatrixType &mattype) const): ditto. + * fMatrix.cc (float rcond): Replace with float rcon everywhere + to avoid shadowed variable warning + (float Matrix::rcond (void) const): New method for + reciprocal condition number calculation. + (float Matrix::rcond (MatrixType &mattype) const): ditto. + * fMatrix.h (float rcond): Replace with float rcon everywhere + to avoid shadowed variable warning + (float Matrix::rcond (void) const): New method for + reciprocal condition number calculation. + (float Matrix::rcond (MatrixType &mattype) const): ditto. + * Array.cc: Fix transpose tests. * CmplxGEBAL.cc (ComplexGEPBALANCE), dbleGEPBAL.cc (GEPBALANCE),
--- a/liboctave/fCMatrix.cc Tue May 13 12:40:23 2008 +0200 +++ b/liboctave/fCMatrix.cc Wed May 14 23:28:41 2008 +0200 @@ -961,45 +961,45 @@ FloatComplexMatrix::inverse (void) const { octave_idx_type info; - float rcond; + float rcon; MatrixType mattype (*this); - return inverse (mattype, info, rcond, 0, 0); + return inverse (mattype, info, rcon, 0, 0); } FloatComplexMatrix FloatComplexMatrix::inverse (octave_idx_type& info) const { - float rcond; + float rcon; MatrixType mattype (*this); - return inverse (mattype, info, rcond, 0, 0); + return inverse (mattype, info, rcon, 0, 0); } FloatComplexMatrix -FloatComplexMatrix::inverse (octave_idx_type& info, float& rcond, int force, +FloatComplexMatrix::inverse (octave_idx_type& info, float& 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); } FloatComplexMatrix FloatComplexMatrix::inverse (MatrixType &mattype) const { octave_idx_type info; - float rcond; - return inverse (mattype, info, rcond, 0, 0); + float rcon; + return inverse (mattype, info, rcon, 0, 0); } FloatComplexMatrix FloatComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info) const { - float rcond; - return inverse (mattype, info, rcond, 0, 0); + float rcon; + return inverse (mattype, info, rcon, 0, 0); } FloatComplexMatrix FloatComplexMatrix::tinverse (MatrixType &mattype, octave_idx_type& info, - float& rcond, int force, int calc_cond) const + float& rcon, int force, int calc_cond) const { FloatComplexMatrix retval; @@ -1023,7 +1023,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) @@ -1037,7 +1037,7 @@ F77_XFCN (ctrcon, CTRCON, (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) @@ -1056,7 +1056,7 @@ FloatComplexMatrix FloatComplexMatrix::finverse (MatrixType &mattype, octave_idx_type& info, - float& rcond, int force, int calc_cond) const + float& rcon, int force, int calc_cond) const { FloatComplexMatrix retval; @@ -1096,7 +1096,7 @@ F77_XFCN (cgetrf, CGETRF, (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) @@ -1108,7 +1108,7 @@ float *prz = rz.fortran_vec (); F77_XFCN (cgecon, CGECON, (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) @@ -1137,7 +1137,7 @@ FloatComplexMatrix FloatComplexMatrix::inverse (MatrixType &mattype, octave_idx_type& info, - float& rcond, int force, int calc_cond) const + float& rcon, int force, int calc_cond) const { int typ = mattype.type (false); FloatComplexMatrix ret; @@ -1146,7 +1146,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 ()) @@ -1155,9 +1155,9 @@ if (info == 0) { if (calc_cond) - rcond = chol.rcond(); + rcon = chol.rcond(); else - rcond = 1.0; + rcon = 1.0; ret = chol.inverse (); } else @@ -1165,9 +1165,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 = FloatComplexMatrix (rows (), columns (), FloatComplex (octave_Float_Inf, 0.)); } @@ -1531,19 +1531,19 @@ FloatComplexMatrix::determinant (void) const { octave_idx_type info; - float rcond; - return determinant (info, rcond, 0); + float rcon; + return determinant (info, rcon, 0); } FloatComplexDET FloatComplexMatrix::determinant (octave_idx_type& info) const { - float rcond; - return determinant (info, rcond, 0); + float rcon; + return determinant (info, rcon, 0); } FloatComplexDET -FloatComplexMatrix::determinant (octave_idx_type& info, float& rcond, int calc_cond) const +FloatComplexMatrix::determinant (octave_idx_type& info, float& rcon, int calc_cond) const { FloatComplexDET retval; @@ -1572,7 +1572,7 @@ F77_XFCN (cgetrf, CGETRF, (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; @@ -1591,7 +1591,7 @@ F77_XFCN (cgecon, CGECON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, tmp_data, nr, anorm, - rcond, pz, prz, info + rcon, pz, prz, info F77_CHAR_ARG_LEN (1))); } @@ -1636,9 +1636,176 @@ return retval; } +float +FloatComplexMatrix::rcond (void) const +{ + MatrixType mattype (*this); + return rcond (mattype); +} + +float +FloatComplexMatrix::rcond (MatrixType &mattype) const +{ + float 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 FloatComplex *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'U'; + char dia = 'N'; + + Array<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (ctrcon, CTRCON, (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 FloatComplex *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'L'; + char dia = 'N'; + + Array<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (ctrcon, CTRCON, (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) + { + float anorm = -1.0; + FloatComplexMatrix atmp = *this; + FloatComplex *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 (cpotrf, CPOTRF, (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<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (cpocon, CPOCON, (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<FloatComplex> z (2 * nc); + FloatComplex *pz = z.fortran_vec (); + Array<float> rz (2 * nc); + float *prz = rz.fortran_vec (); + + F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info)); + + if (info != 0) + { + rcon = 0.0; + mattype.mark_as_rectangular (); + } + else + { + char job = '1'; + F77_XFCN (cgecon, CGECON, (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; +} + FloatComplexMatrix FloatComplexMatrix::utsolve (MatrixType &mattype, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { @@ -1660,7 +1827,7 @@ typ == MatrixType::Upper) { octave_idx_type b_nc = b.cols (); - rcond = 1.; + rcon = 1.; info = 0; if (typ == MatrixType::Permuted_Upper) @@ -1686,7 +1853,7 @@ F77_XFCN (ctrcon, CTRCON, (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) @@ -1695,18 +1862,18 @@ if (info != 0) info = -2; - volatile float rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile float 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); } } @@ -1739,7 +1906,7 @@ FloatComplexMatrix FloatComplexMatrix::ltsolve (MatrixType &mattype, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { @@ -1761,7 +1928,7 @@ typ == MatrixType::Lower) { octave_idx_type b_nc = b.cols (); - rcond = 1.; + rcon = 1.; info = 0; if (typ == MatrixType::Permuted_Lower) @@ -1787,7 +1954,7 @@ F77_XFCN (ctrcon, CTRCON, (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) @@ -1796,18 +1963,18 @@ if (info != 0) info = -2; - volatile float rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile float 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); } } @@ -1840,7 +2007,7 @@ FloatComplexMatrix FloatComplexMatrix::fsolve (MatrixType &mattype, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { @@ -1875,7 +2042,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; @@ -1894,24 +2061,24 @@ F77_XFCN (cpocon, CPOCON, (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 float rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile float 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); } } @@ -1957,13 +2124,13 @@ F77_XFCN (cgetrf, CGETRF, (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"); @@ -1979,24 +2146,24 @@ char job = '1'; F77_XFCN (cgecon, CGECON, (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 float rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile float 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); } } @@ -2026,60 +2193,60 @@ FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b) const { octave_idx_type info; - float rcond; - return solve (typ, b, info, rcond, 0); + float rcon; + return solve (typ, b, info, rcon, 0); } FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info) const { - float rcond; - return solve (typ, b, info, rcond, 0); + float rcon; + return solve (typ, b, info, rcon, 0); } FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcond) const + float& rcon) const { - return solve (typ, b, info, rcond, 0); + return solve (typ, b, info, rcon, 0); } FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler, + float& rcon, solve_singularity_handler sing_handler, bool singular_fallback) const { FloatComplexMatrix tmp (b); - return solve (typ, tmp, info, rcond, sing_handler, singular_fallback); + return solve (typ, tmp, info, rcon, sing_handler, singular_fallback); } FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b) const { octave_idx_type info; - float rcond; - return solve (typ, b, info, rcond, 0); + float rcon; + return solve (typ, b, info, rcon, 0); } FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info) const { - float rcond; - return solve (typ, b, info, rcond, 0); + float rcon; + return solve (typ, b, info, rcon, 0); } FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond) const + octave_idx_type& info, float& rcon) const { - return solve (typ, b, info, rcond, 0); + return solve (typ, b, info, rcon, 0); } FloatComplexMatrix FloatComplexMatrix::solve (MatrixType &mattype, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, bool singular_fallback) const { @@ -2091,11 +2258,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"); @@ -2106,7 +2273,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; @@ -2116,183 +2283,183 @@ FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b) const { octave_idx_type info; - float rcond; - return solve (typ, FloatComplexColumnVector (b), info, rcond, 0); + float rcon; + return solve (typ, FloatComplexColumnVector (b), info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info) const { - float rcond; - return solve (typ, FloatComplexColumnVector (b), info, rcond, 0); + float rcon; + return solve (typ, FloatComplexColumnVector (b), info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b, - octave_idx_type& info, float& rcond) const + octave_idx_type& info, float& rcon) const { - return solve (typ, FloatComplexColumnVector (b), info, rcond, 0); + return solve (typ, FloatComplexColumnVector (b), info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatColumnVector& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const { - return solve (typ, FloatComplexColumnVector (b), info, rcond, sing_handler); + return solve (typ, FloatComplexColumnVector (b), info, rcon, sing_handler); } FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b) const { octave_idx_type info; - float rcond; - return solve (typ, b, info, rcond, 0); + float rcon; + return solve (typ, b, info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, octave_idx_type& info) const { - float rcond; - return solve (typ, b, info, rcond, 0); + float rcon; + return solve (typ, b, info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info, float& rcond) const + octave_idx_type& info, float& rcon) const { - return solve (typ, b, info, rcond, 0); + return solve (typ, b, info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const { FloatComplexMatrix 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)); } FloatComplexMatrix FloatComplexMatrix::solve (const FloatMatrix& b) const { octave_idx_type info; - float rcond; - return solve (b, info, rcond, 0); + float rcon; + return solve (b, info, rcon, 0); } FloatComplexMatrix FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info) const { - float rcond; - return solve (b, info, rcond, 0); + float rcon; + return solve (b, info, rcon, 0); } FloatComplexMatrix -FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info, float& rcond) const +FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info, float& rcon) const { - return solve (b, info, rcond, 0); + return solve (b, info, rcon, 0); } FloatComplexMatrix -FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info, float& rcond, +FloatComplexMatrix::solve (const FloatMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const { FloatComplexMatrix tmp (b); - return solve (tmp, info, rcond, sing_handler); + return solve (tmp, info, rcon, sing_handler); } FloatComplexMatrix FloatComplexMatrix::solve (const FloatComplexMatrix& b) const { octave_idx_type info; - float rcond; - return solve (b, info, rcond, 0); + float rcon; + return solve (b, info, rcon, 0); } FloatComplexMatrix FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info) const { - float rcond; - return solve (b, info, rcond, 0); + float rcon; + return solve (b, info, rcon, 0); } FloatComplexMatrix -FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcond) const +FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon) const { - return solve (b, info, rcond, 0); + return solve (b, info, rcon, 0); } FloatComplexMatrix -FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcond, +FloatComplexMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& 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); } FloatComplexColumnVector FloatComplexMatrix::solve (const FloatColumnVector& b) const { octave_idx_type info; - float rcond; - return solve (FloatComplexColumnVector (b), info, rcond, 0); + float rcon; + return solve (FloatComplexColumnVector (b), info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (const FloatColumnVector& b, octave_idx_type& info) const { - float rcond; - return solve (FloatComplexColumnVector (b), info, rcond, 0); + float rcon; + return solve (FloatComplexColumnVector (b), info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, - float& rcond) const + float& rcon) const { - return solve (FloatComplexColumnVector (b), info, rcond, 0); + return solve (FloatComplexColumnVector (b), info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, - float& rcond, + float& rcon, solve_singularity_handler sing_handler) const { - return solve (FloatComplexColumnVector (b), info, rcond, sing_handler); + return solve (FloatComplexColumnVector (b), info, rcon, sing_handler); } FloatComplexColumnVector FloatComplexMatrix::solve (const FloatComplexColumnVector& b) const { octave_idx_type info; - float rcond; - return solve (b, info, rcond, 0); + float rcon; + return solve (b, info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info) const { - float rcond; - return solve (b, info, rcond, 0); + float rcon; + return solve (b, info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, - float& rcond) const + float& rcon) const { - return solve (b, info, rcond, 0); + return solve (b, info, rcon, 0); } FloatComplexColumnVector FloatComplexMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, - float& rcond, + float& 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); } FloatComplexMatrix @@ -2300,31 +2467,31 @@ { octave_idx_type info; octave_idx_type rank; - float rcond; - return lssolve (FloatComplexMatrix (b), info, rank, rcond); + float rcon; + return lssolve (FloatComplexMatrix (b), info, rank, rcon); } FloatComplexMatrix FloatComplexMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info) const { octave_idx_type rank; - float rcond; - return lssolve (FloatComplexMatrix (b), info, rank, rcond); + float rcon; + return lssolve (FloatComplexMatrix (b), info, rank, rcon); } FloatComplexMatrix FloatComplexMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info, octave_idx_type& rank) const { - float rcond; - return lssolve (FloatComplexMatrix (b), info, rank, rcond); + float rcon; + return lssolve (FloatComplexMatrix (b), info, rank, rcon); } FloatComplexMatrix FloatComplexMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const + octave_idx_type& rank, float& rcon) const { - return lssolve (FloatComplexMatrix (b), info, rank, rcond); + return lssolve (FloatComplexMatrix (b), info, rank, rcon); } FloatComplexMatrix @@ -2332,29 +2499,29 @@ { octave_idx_type info; octave_idx_type rank; - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatComplexMatrix FloatComplexMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info) const { octave_idx_type rank; - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatComplexMatrix FloatComplexMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info, octave_idx_type& rank) const { - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatComplexMatrix FloatComplexMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const + octave_idx_type& rank, float& rcon) const { FloatComplexMatrix retval; @@ -2372,7 +2539,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) { @@ -2439,7 +2606,7 @@ octave_idx_type* piwork = iwork.fortran_vec (); F77_XFCN (cgelsd, CGELSD, (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 @@ -2476,19 +2643,19 @@ work.resize (lwork); F77_XFCN (cgelsd, CGELSD, (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); } @@ -2501,31 +2668,31 @@ { octave_idx_type info; octave_idx_type rank; - float rcond; - return lssolve (FloatComplexColumnVector (b), info, rank, rcond); + float rcon; + return lssolve (FloatComplexColumnVector (b), info, rank, rcon); } FloatComplexColumnVector FloatComplexMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info) const { octave_idx_type rank; - float rcond; - return lssolve (FloatComplexColumnVector (b), info, rank, rcond); + float rcon; + return lssolve (FloatComplexColumnVector (b), info, rank, rcon); } FloatComplexColumnVector FloatComplexMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const { - float rcond; - return lssolve (FloatComplexColumnVector (b), info, rank, rcond); + float rcon; + return lssolve (FloatComplexColumnVector (b), info, rank, rcon); } FloatComplexColumnVector FloatComplexMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const + octave_idx_type& rank, float& rcon) const { - return lssolve (FloatComplexColumnVector (b), info, rank, rcond); + return lssolve (FloatComplexColumnVector (b), info, rank, rcon); } FloatComplexColumnVector @@ -2533,30 +2700,30 @@ { octave_idx_type info; octave_idx_type rank; - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatComplexColumnVector FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info) const { octave_idx_type rank; - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatComplexColumnVector FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const { - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatComplexColumnVector FloatComplexMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const + octave_idx_type& rank, float& rcon) const { FloatComplexColumnVector retval; @@ -2574,7 +2741,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) { @@ -2633,7 +2800,7 @@ octave_idx_type* piwork = iwork.fortran_vec (); F77_XFCN (cgelsd, CGELSD, (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))); @@ -2642,7 +2809,7 @@ iwork.resize (iwork(0)); F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcond, rank, + maxmn, ps, rcon, rank, work.fortran_vec (), lwork, prwork, piwork, info)); @@ -2651,12 +2818,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); } @@ -2680,11 +2847,11 @@ }; static void -solve_singularity_warning (float rcond) +solve_singularity_warning (float rcon) { (*current_liboctave_warning_handler) ("singular matrix encountered in expm calculation, rcond = %g", - rcond); + rcon); } FloatComplexMatrix @@ -2811,8 +2978,8 @@ // Compute pade approximation = inverse (dpp) * npp. - float rcond; - retval = dpp.solve (npp, info, rcond, solve_singularity_warning); + float rcon; + retval = dpp.solve (npp, info, rcon, solve_singularity_warning); if (info < 0) return retval;
--- a/liboctave/fCMatrix.h Tue May 13 12:40:23 2008 +0200 +++ b/liboctave/fCMatrix.h Wed May 14 23:28:41 2008 +0200 @@ -38,7 +38,7 @@ { public: - typedef void (*solve_singularity_handler) (float rcond); + typedef void (*solve_singularity_handler) (float rcon); FloatComplexMatrix (void) : MArray2<FloatComplex> () { } @@ -147,21 +147,21 @@ private: FloatComplexMatrix tinverse (MatrixType &mattype, octave_idx_type& info, - float& rcond, int force, int calc_cond) const; + float& rcon, int force, int calc_cond) const; FloatComplexMatrix finverse (MatrixType &mattype, octave_idx_type& info, - float& rcond, int force, int calc_cond) const; + float& rcon, int force, int calc_cond) const; public: FloatComplexMatrix inverse (void) const; FloatComplexMatrix inverse (octave_idx_type& info) const; - FloatComplexMatrix inverse (octave_idx_type& info, float& rcond, int force = 0, + FloatComplexMatrix inverse (octave_idx_type& info, float& rcon, int force = 0, int calc_cond = 1) const; FloatComplexMatrix inverse (MatrixType &mattype) const; FloatComplexMatrix inverse (MatrixType &mattype, octave_idx_type& info) const; FloatComplexMatrix inverse (MatrixType &mattype, octave_idx_type& info, - float& rcond, int force = 0, + float& rcon, int force = 0, int calc_cond = 1) const; FloatComplexMatrix pseudo_inverse (float tol = 0.0) const; @@ -174,24 +174,27 @@ FloatComplexDET determinant (void) const; FloatComplexDET determinant (octave_idx_type& info) const; - FloatComplexDET determinant (octave_idx_type& info, float& rcond, int calc_cond = 1) const; + FloatComplexDET determinant (octave_idx_type& info, float& rcon, int calc_cond = 1) const; + + float rcond (void) const; + float rcond (MatrixType &mattype) const; private: // Upper triangular matrix solvers FloatComplexMatrix utsolve (MatrixType &typ, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; // Lower triangular matrix solvers FloatComplexMatrix ltsolve (MatrixType &typ, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; // Full matrix solvers (umfpack/cholesky) FloatComplexMatrix fsolve (MatrixType &typ, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; @@ -201,18 +204,18 @@ FloatComplexMatrix solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info) const; FloatComplexMatrix solve (MatrixType &typ, const FloatMatrix& b, - octave_idx_type& info, float& rcond) const; + octave_idx_type& info, float& rcon) const; FloatComplexMatrix solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler, + float& rcon, solve_singularity_handler sing_handler, bool singular_fallback = true) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond) const; + octave_idx_type& info, float& rcon) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, bool singular_fallback = true) const; @@ -220,9 +223,9 @@ FloatComplexColumnVector solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatColumnVector& b, - octave_idx_type& info, float& rcond) const; + octave_idx_type& info, float& rcon) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatColumnVector& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const; FloatComplexColumnVector solve (MatrixType &typ, @@ -230,37 +233,37 @@ FloatComplexColumnVector solve (MatrixType &typ, const FloatComplexColumnVector& b, octave_idx_type& info) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info, float& rcond) const; + octave_idx_type& info, float& rcon) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const; // Generic interface to solver with probing of type FloatComplexMatrix solve (const FloatMatrix& b) const; FloatComplexMatrix solve (const FloatMatrix& b, octave_idx_type& info) const; - FloatComplexMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcond) const; - FloatComplexMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcond, + FloatComplexMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcon) const; + FloatComplexMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const; FloatComplexMatrix solve (const FloatComplexMatrix& b) const; FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info) const; - FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcond) const; - FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcond, + FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon) const; + FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const; FloatComplexColumnVector solve (const FloatColumnVector& b) const; FloatComplexColumnVector solve (const FloatColumnVector& b, octave_idx_type& info) const; FloatComplexColumnVector solve (const FloatColumnVector& b, octave_idx_type& info, - float& rcond) const; - FloatComplexColumnVector solve (const FloatColumnVector& b, octave_idx_type& info, float& rcond, + float& rcon) const; + FloatComplexColumnVector solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b, octave_idx_type& info) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b, octave_idx_type& info, - float& rcond) const; + float& rcon) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b, octave_idx_type& info, - float& rcond, + float& rcon, solve_singularity_handler sing_handler) const; FloatComplexMatrix lssolve (const FloatMatrix& b) const; @@ -268,14 +271,14 @@ FloatComplexMatrix lssolve (const FloatMatrix& b, octave_idx_type& info, octave_idx_type& rank) const; FloatComplexMatrix lssolve (const FloatMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const; + octave_idx_type& rank, float& rcon) const; FloatComplexMatrix lssolve (const FloatComplexMatrix& b) const; FloatComplexMatrix lssolve (const FloatComplexMatrix& b, octave_idx_type& info) const; FloatComplexMatrix lssolve (const FloatComplexMatrix& b, octave_idx_type& info, octave_idx_type& rank) const; FloatComplexMatrix lssolve (const FloatComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const; + octave_idx_type& rank, float& rcon) const; FloatComplexColumnVector lssolve (const FloatColumnVector& b) const; FloatComplexColumnVector lssolve (const FloatColumnVector& b, @@ -283,7 +286,7 @@ FloatComplexColumnVector lssolve (const FloatColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const; FloatComplexColumnVector lssolve (const FloatColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const; + octave_idx_type& rank, float& rcon) const; FloatComplexColumnVector lssolve (const FloatComplexColumnVector& b) const; FloatComplexColumnVector lssolve (const FloatComplexColumnVector& b, @@ -293,7 +296,7 @@ octave_idx_type& rank) const; FloatComplexColumnVector lssolve (const FloatComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const; + octave_idx_type& rank, float& rcon) const; FloatComplexMatrix expm (void) const;
--- a/liboctave/fMatrix.cc Tue May 13 12:40:23 2008 +0200 +++ b/liboctave/fMatrix.cc Wed May 14 23:28:41 2008 +0200 @@ -650,44 +650,44 @@ FloatMatrix::inverse (void) const { octave_idx_type info; - float rcond; + float rcon; MatrixType mattype (*this); - return inverse (mattype, info, rcond, 0, 0); + return inverse (mattype, info, rcon, 0, 0); } FloatMatrix FloatMatrix::inverse (octave_idx_type& info) const { - float rcond; + float rcon; MatrixType mattype (*this); - return inverse (mattype, info, rcond, 0, 0); + return inverse (mattype, info, rcon, 0, 0); } FloatMatrix -FloatMatrix::inverse (octave_idx_type& info, float& rcond, int force, +FloatMatrix::inverse (octave_idx_type& info, float& 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); } FloatMatrix FloatMatrix::inverse (MatrixType& mattype) const { octave_idx_type info; - float rcond; - return inverse (mattype, info, rcond, 0, 0); + float rcon; + return inverse (mattype, info, rcon, 0, 0); } FloatMatrix FloatMatrix::inverse (MatrixType &mattype, octave_idx_type& info) const { - float rcond; - return inverse (mattype, info, rcond, 0, 0); + float rcon; + return inverse (mattype, info, rcon, 0, 0); } FloatMatrix -FloatMatrix::tinverse (MatrixType &mattype, octave_idx_type& info, float& rcond, +FloatMatrix::tinverse (MatrixType &mattype, octave_idx_type& info, float& rcon, int force, int calc_cond) const { FloatMatrix retval; @@ -712,7 +712,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) @@ -726,7 +726,7 @@ F77_XFCN (strcon, STRCON, (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) @@ -745,7 +745,7 @@ FloatMatrix -FloatMatrix::finverse (MatrixType &mattype, octave_idx_type& info, float& rcond, +FloatMatrix::finverse (MatrixType &mattype, octave_idx_type& info, float& rcon, int force, int calc_cond) const { FloatMatrix retval; @@ -785,7 +785,7 @@ F77_XFCN (sgetrf, SGETRF, (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) @@ -798,7 +798,7 @@ octave_idx_type *piz = iz.fortran_vec (); F77_XFCN (sgecon, SGECON, (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) @@ -826,7 +826,7 @@ } FloatMatrix -FloatMatrix::inverse (MatrixType &mattype, octave_idx_type& info, float& rcond, +FloatMatrix::inverse (MatrixType &mattype, octave_idx_type& info, float& rcon, int force, int calc_cond) const { int typ = mattype.type (false); @@ -836,7 +836,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 ()) @@ -845,9 +845,9 @@ if (info == 0) { if (calc_cond) - rcond = chol.rcond (); + rcon = chol.rcond (); else - rcond = 1.0; + rcon = 1.0; ret = chol.inverse (); } else @@ -855,9 +855,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 = FloatMatrix (rows (), columns (), octave_Float_Inf); } @@ -1215,19 +1215,19 @@ FloatMatrix::determinant (void) const { octave_idx_type info; - float rcond; - return determinant (info, rcond, 0); + float rcon; + return determinant (info, rcon, 0); } FloatDET FloatMatrix::determinant (octave_idx_type& info) const { - float rcond; - return determinant (info, rcond, 0); + float rcon; + return determinant (info, rcon, 0); } FloatDET -FloatMatrix::determinant (octave_idx_type& info, float& rcond, int calc_cond) const +FloatMatrix::determinant (octave_idx_type& info, float& rcon, int calc_cond) const { FloatDET retval; @@ -1256,7 +1256,7 @@ F77_XFCN (sgetrf, SGETRF, (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; @@ -1275,7 +1275,7 @@ F77_XFCN (sgecon, SGECON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, tmp_data, nr, anorm, - rcond, pz, piz, info + rcon, pz, piz, info F77_CHAR_ARG_LEN (1))); } @@ -1320,9 +1320,174 @@ return retval; } +float +FloatMatrix::rcond (void) const +{ + MatrixType mattype (*this); + return rcond (mattype); +} + +float +FloatMatrix::rcond (MatrixType &mattype) const +{ + float 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 float *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'U'; + char dia = 'N'; + + Array<float> z (3 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (strcon, STRCON, (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 float *tmp_data = fortran_vec (); + octave_idx_type info = 0; + char norm = '1'; + char uplo = 'L'; + char dia = 'N'; + + Array<float> z (3 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (strcon, STRCON, (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) + { + float anorm = -1.0; + FloatMatrix atmp = *this; + float *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 (spotrf, SPOTRF, (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<float> z (3 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (spocon, SPOCON, (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<float> z (4 * nc); + float *pz = z.fortran_vec (); + Array<octave_idx_type> iz (nc); + octave_idx_type *piz = iz.fortran_vec (); + + F77_XFCN (sgetrf, SGETRF, (nr, nr, tmp_data, nr, pipvt, info)); + + if (info != 0) + { + rcon = 0.0; + mattype.mark_as_rectangular (); + } + else + { + char job = '1'; + F77_XFCN (sgecon, SGECON, (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; +} + FloatMatrix FloatMatrix::utsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler, + float& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { FloatMatrix retval; @@ -1343,7 +1508,7 @@ typ == MatrixType::Upper) { octave_idx_type b_nc = b.cols (); - rcond = 1.; + rcon = 1.; info = 0; if (typ == MatrixType::Permuted_Upper) @@ -1369,7 +1534,7 @@ F77_XFCN (strcon, STRCON, (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) @@ -1378,18 +1543,18 @@ if (info != 0) info = -2; - volatile float rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile float 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); } } @@ -1422,7 +1587,7 @@ FloatMatrix FloatMatrix::ltsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler, + float& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { FloatMatrix retval; @@ -1443,7 +1608,7 @@ typ == MatrixType::Lower) { octave_idx_type b_nc = b.cols (); - rcond = 1.; + rcon = 1.; info = 0; if (typ == MatrixType::Permuted_Lower) @@ -1469,7 +1634,7 @@ F77_XFCN (strcon, STRCON, (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) @@ -1478,18 +1643,18 @@ if (info != 0) info = -2; - volatile float rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile float 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); } } @@ -1522,7 +1687,7 @@ FloatMatrix FloatMatrix::fsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler, + float& rcon, solve_singularity_handler sing_handler, bool calc_cond) const { FloatMatrix retval; @@ -1555,7 +1720,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; @@ -1574,24 +1739,24 @@ F77_XFCN (spocon, SPOCON, (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 float rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile float 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); } } @@ -1635,13 +1800,13 @@ F77_XFCN (sgetrf, SGETRF, (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"); @@ -1657,24 +1822,24 @@ char job = '1'; F77_XFCN (sgecon, SGECON, (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 float rcond_plus_one = rcond + 1.0; - - if (rcond_plus_one == 1.0 || xisnan (rcond)) + volatile float 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); } } @@ -1706,20 +1871,20 @@ FloatMatrix::solve (MatrixType &typ, const FloatMatrix& b) const { octave_idx_type info; - float rcond; - return solve (typ, b, info, rcond, 0); + float rcon; + return solve (typ, b, info, rcon, 0); } FloatMatrix FloatMatrix::solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcond) const + float& rcon) const { - return solve (typ, b, info, rcond, 0); + return solve (typ, b, info, rcon, 0); } FloatMatrix FloatMatrix::solve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler, + float& rcon, solve_singularity_handler sing_handler, bool singular_fallback) const { FloatMatrix retval; @@ -1730,11 +1895,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"); @@ -1745,7 +1910,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; @@ -1768,49 +1933,49 @@ FloatComplexMatrix FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info, - float& rcond) const + float& rcon) const { FloatComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcond); + return tmp.solve (typ, b, info, rcon); } FloatComplexMatrix FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler, + float& rcon, solve_singularity_handler sing_handler, bool singular_fallback) const { FloatComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcond, sing_handler, singular_fallback); + return tmp.solve (typ, b, info, rcon, sing_handler, singular_fallback); } FloatColumnVector FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b) const { - octave_idx_type info; float rcond; - return solve (typ, b, info, rcond); + octave_idx_type info; float rcon; + return solve (typ, b, info, rcon); } FloatColumnVector FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info) const { - float rcond; - return solve (typ, b, info, rcond); + float rcon; + return solve (typ, b, info, rcon); } FloatColumnVector FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info, - float& rcond) const + float& rcon) const { - return solve (typ, b, info, rcond, 0); + return solve (typ, b, info, rcon, 0); } FloatColumnVector FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler) const + float& rcon, solve_singularity_handler sing_handler) const { FloatMatrix 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)); } FloatComplexColumnVector @@ -1830,48 +1995,48 @@ FloatComplexColumnVector FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info, float& rcond) const + octave_idx_type& info, float& rcon) const { FloatComplexMatrix tmp (*this); - return tmp.solve (typ, b, info, rcond); + return tmp.solve (typ, b, info, rcon); } FloatComplexColumnVector FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const { FloatComplexMatrix tmp (*this); - return tmp.solve(typ, b, info, rcond, sing_handler); + return tmp.solve(typ, b, info, rcon, sing_handler); } FloatMatrix FloatMatrix::solve (const FloatMatrix& b) const { octave_idx_type info; - float rcond; - return solve (b, info, rcond, 0); + float rcon; + return solve (b, info, rcon, 0); } FloatMatrix FloatMatrix::solve (const FloatMatrix& b, octave_idx_type& info) const { - float rcond; - return solve (b, info, rcond, 0); + float rcon; + return solve (b, info, rcon, 0); } FloatMatrix -FloatMatrix::solve (const FloatMatrix& b, octave_idx_type& info, float& rcond) const +FloatMatrix::solve (const FloatMatrix& b, octave_idx_type& info, float& rcon) const { - return solve (b, info, rcond, 0); + return solve (b, info, rcon, 0); } FloatMatrix FloatMatrix::solve (const FloatMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler) const + float& 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); } FloatComplexMatrix @@ -1889,46 +2054,46 @@ } FloatComplexMatrix -FloatMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcond) const +FloatMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon) const { FloatComplexMatrix tmp (*this); - return tmp.solve (b, info, rcond); + return tmp.solve (b, info, rcon); } FloatComplexMatrix -FloatMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcond, +FloatMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const { FloatComplexMatrix tmp (*this); - return tmp.solve (b, info, rcond, sing_handler); + return tmp.solve (b, info, rcon, sing_handler); } FloatColumnVector FloatMatrix::solve (const FloatColumnVector& b) const { - octave_idx_type info; float rcond; - return solve (b, info, rcond); + octave_idx_type info; float rcon; + return solve (b, info, rcon); } FloatColumnVector FloatMatrix::solve (const FloatColumnVector& b, octave_idx_type& info) const { - float rcond; - return solve (b, info, rcond); + float rcon; + return solve (b, info, rcon); } FloatColumnVector -FloatMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, float& rcond) const +FloatMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon) const { - return solve (b, info, rcond, 0); + return solve (b, info, rcon, 0); } FloatColumnVector -FloatMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, float& rcond, +FloatMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, float& 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); } FloatComplexColumnVector @@ -1946,18 +2111,18 @@ } FloatComplexColumnVector -FloatMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, float& rcond) const +FloatMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, float& rcon) const { FloatComplexMatrix tmp (*this); - return tmp.solve (b, info, rcond); + return tmp.solve (b, info, rcon); } FloatComplexColumnVector -FloatMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, float& rcond, +FloatMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const { FloatComplexMatrix tmp (*this); - return tmp.solve (b, info, rcond, sing_handler); + return tmp.solve (b, info, rcon, sing_handler); } FloatMatrix @@ -1965,29 +2130,29 @@ { octave_idx_type info; octave_idx_type rank; - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatMatrix FloatMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info) const { octave_idx_type rank; - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatMatrix FloatMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info, octave_idx_type& rank) const { - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatMatrix FloatMatrix::lssolve (const FloatMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float &rcond) const + octave_idx_type& rank, float &rcon) const { FloatMatrix retval; @@ -2005,7 +2170,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 = FloatMatrix (maxmn, nrhs, 0.0); @@ -2063,7 +2228,7 @@ octave_idx_type* piwork = iwork.fortran_vec (); F77_XFCN (sgelsd, SGELSD, (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 @@ -2107,7 +2272,7 @@ work.resize (lwork); F77_XFCN (sgelsd, SGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcond, rank, + maxmn, ps, rcon, rank, work.fortran_vec (), lwork, piwork, info)); @@ -2115,9 +2280,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); } @@ -2131,8 +2296,8 @@ FloatComplexMatrix tmp (*this); octave_idx_type info; octave_idx_type rank; - float rcond; - return tmp.lssolve (b, info, rank, rcond); + float rcon; + return tmp.lssolve (b, info, rank, rcon); } FloatComplexMatrix @@ -2140,8 +2305,8 @@ { FloatComplexMatrix tmp (*this); octave_idx_type rank; - float rcond; - return tmp.lssolve (b, info, rank, rcond); + float rcon; + return tmp.lssolve (b, info, rank, rcon); } FloatComplexMatrix @@ -2149,16 +2314,16 @@ octave_idx_type& rank) const { FloatComplexMatrix tmp (*this); - float rcond; - return tmp.lssolve (b, info, rank, rcond); + float rcon; + return tmp.lssolve (b, info, rank, rcon); } FloatComplexMatrix FloatMatrix::lssolve (const FloatComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const + octave_idx_type& rank, float& rcon) const { FloatComplexMatrix tmp (*this); - return tmp.lssolve (b, info, rank, rcond); + return tmp.lssolve (b, info, rank, rcon); } FloatColumnVector @@ -2166,29 +2331,29 @@ { octave_idx_type info; octave_idx_type rank; - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatColumnVector FloatMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info) const { octave_idx_type rank; - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatColumnVector FloatMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const { - float rcond; - return lssolve (b, info, rank, rcond); + float rcon; + return lssolve (b, info, rank, rcon); } FloatColumnVector FloatMatrix::lssolve (const FloatColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float &rcond) const + octave_idx_type& rank, float &rcon) const { FloatColumnVector retval; @@ -2206,7 +2371,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) { @@ -2257,14 +2422,14 @@ octave_idx_type* piwork = iwork.fortran_vec (); F77_XFCN (sgelsd, SGELSD, (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 (sgelsd, SGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcond, rank, + maxmn, ps, rcon, rank, work.fortran_vec (), lwork, piwork, info)); @@ -2274,9 +2439,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); @@ -2291,8 +2456,8 @@ FloatComplexMatrix tmp (*this); octave_idx_type info; octave_idx_type rank; - float rcond; - return tmp.lssolve (b, info, rank, rcond); + float rcon; + return tmp.lssolve (b, info, rank, rcon); } FloatComplexColumnVector @@ -2300,8 +2465,8 @@ { FloatComplexMatrix tmp (*this); octave_idx_type rank; - float rcond; - return tmp.lssolve (b, info, rank, rcond); + float rcon; + return tmp.lssolve (b, info, rank, rcon); } FloatComplexColumnVector @@ -2309,16 +2474,16 @@ octave_idx_type& rank) const { FloatComplexMatrix tmp (*this); - float rcond; - return tmp.lssolve (b, info, rank, rcond); + float rcon; + return tmp.lssolve (b, info, rank, rcon); } FloatComplexColumnVector FloatMatrix::lssolve (const FloatComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float &rcond) const + octave_idx_type& rank, float &rcon) const { FloatComplexMatrix tmp (*this); - return tmp.lssolve (b, info, rank, rcond); + return tmp.lssolve (b, info, rank, rcon); } // Constants for matrix exponential calculation. @@ -2336,11 +2501,11 @@ }; static void -solve_singularity_warning (float rcond) +solve_singularity_warning (float rcon) { (*current_liboctave_warning_handler) ("singular matrix encountered in expm calculation, rcond = %g", - rcond); + rcon); } FloatMatrix @@ -2465,8 +2630,8 @@ // Compute pade approximation = inverse (dpp) * npp. - float rcond; - retval = dpp.solve (npp, info, rcond, solve_singularity_warning); + float rcon; + retval = dpp.solve (npp, info, rcon, solve_singularity_warning); if (info < 0) return retval;
--- a/liboctave/fMatrix.h Tue May 13 12:40:23 2008 +0200 +++ b/liboctave/fMatrix.h Wed May 14 23:28:41 2008 +0200 @@ -37,7 +37,7 @@ { public: - typedef void (*solve_singularity_handler) (float rcond); + typedef void (*solve_singularity_handler) (float rcon); FloatMatrix (void) : MArray2<float> () { } @@ -116,21 +116,21 @@ FloatColumnVector column (octave_idx_type i) const; private: - FloatMatrix tinverse (MatrixType &mattype, octave_idx_type& info, float& rcond, + FloatMatrix tinverse (MatrixType &mattype, octave_idx_type& info, float& rcon, int force, int calc_cond) const; - FloatMatrix finverse (MatrixType &mattype, octave_idx_type& info, float& rcond, + FloatMatrix finverse (MatrixType &mattype, octave_idx_type& info, float& rcon, int force, int calc_cond) const; public: FloatMatrix inverse (void) const; FloatMatrix inverse (octave_idx_type& info) const; - FloatMatrix inverse (octave_idx_type& info, float& rcond, int force = 0, + FloatMatrix inverse (octave_idx_type& info, float& rcon, int force = 0, int calc_cond = 1) const; FloatMatrix inverse (MatrixType &mattype) const; FloatMatrix inverse (MatrixType &mattype, octave_idx_type& info) const; - FloatMatrix inverse (MatrixType &mattype, octave_idx_type& info, float& rcond, + FloatMatrix inverse (MatrixType &mattype, octave_idx_type& info, float& rcon, int force = 0, int calc_cond = 1) const; FloatMatrix pseudo_inverse (float tol = 0.0) const; @@ -143,22 +143,25 @@ FloatDET determinant (void) const; FloatDET determinant (octave_idx_type& info) const; - FloatDET determinant (octave_idx_type& info, float& rcond, int calc_cond = 1) const; + FloatDET determinant (octave_idx_type& info, float& rcon, int calc_cond = 1) const; + + float rcond (void) const; + float rcond (MatrixType &mattype) const; private: // Upper triangular matrix solvers FloatMatrix utsolve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler, + float& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; // Lower triangular matrix solvers FloatMatrix ltsolve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler, + float& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; // Full matrix solvers (lu/cholesky) FloatMatrix fsolve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler, + float& rcon, solve_singularity_handler sing_handler, bool calc_cond = false) const; public: @@ -166,18 +169,18 @@ FloatMatrix solve (MatrixType &typ, const FloatMatrix& b) const; FloatMatrix solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info) const; FloatMatrix solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcond) const; + float& rcon) const; FloatMatrix solve (MatrixType &typ, const FloatMatrix& b, octave_idx_type& info, - float& rcond, solve_singularity_handler sing_handler, + float& rcon, solve_singularity_handler sing_handler, bool singular_fallback = true) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond) const; + octave_idx_type& info, float& rcon) const; FloatComplexMatrix solve (MatrixType &typ, const FloatComplexMatrix& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler, bool singular_fallback = true) const; @@ -185,9 +188,9 @@ FloatColumnVector solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info) const; FloatColumnVector solve (MatrixType &typ, const FloatColumnVector& b, - octave_idx_type& info, float& rcond) const; + octave_idx_type& info, float& rcon) const; FloatColumnVector solve (MatrixType &typ, const FloatColumnVector& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const; FloatComplexColumnVector solve (MatrixType &typ, @@ -195,36 +198,36 @@ FloatComplexColumnVector solve (MatrixType &typ, const FloatComplexColumnVector& b, octave_idx_type& info) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info, float& rcond) const; + octave_idx_type& info, float& rcon) const; FloatComplexColumnVector solve (MatrixType &typ, const FloatComplexColumnVector& b, - octave_idx_type& info, float& rcond, + octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const; // Generic interface to solver with probing of type FloatMatrix solve (const FloatMatrix& b) const; FloatMatrix solve (const FloatMatrix& b, octave_idx_type& info) const; - FloatMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcond) const; - FloatMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcond, + FloatMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcon) const; + FloatMatrix solve (const FloatMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const; FloatComplexMatrix solve (const FloatComplexMatrix& b) const; FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info) const; - FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcond) const; - FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcond, + FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon) const; + FloatComplexMatrix solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const; FloatColumnVector solve (const FloatColumnVector& b) const; FloatColumnVector solve (const FloatColumnVector& b, octave_idx_type& info) const; - FloatColumnVector solve (const FloatColumnVector& b, octave_idx_type& info, float& rcond) const; - FloatColumnVector solve (const FloatColumnVector& b, octave_idx_type& info, float& rcond, + FloatColumnVector solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon) const; + FloatColumnVector solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon, solve_singularity_handler sing_handler) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b, octave_idx_type& info) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b, octave_idx_type& info, - float& rcond) const; + float& rcon) const; FloatComplexColumnVector solve (const FloatComplexColumnVector& b, octave_idx_type& info, - float& rcond, + float& rcon, solve_singularity_handler sing_handler) const; // Singular solvers @@ -233,21 +236,21 @@ FloatMatrix lssolve (const FloatMatrix& b, octave_idx_type& info, octave_idx_type& rank) const; FloatMatrix lssolve (const FloatMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const; + octave_idx_type& rank, float& rcon) const; FloatComplexMatrix lssolve (const FloatComplexMatrix& b) const; FloatComplexMatrix lssolve (const FloatComplexMatrix& b, octave_idx_type& info) const; FloatComplexMatrix lssolve (const FloatComplexMatrix& b, octave_idx_type& info, octave_idx_type& rank) const; FloatComplexMatrix lssolve (const FloatComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank, float &rcond) const; + octave_idx_type& rank, float &rcon) const; FloatColumnVector lssolve (const FloatColumnVector& b) const; FloatColumnVector lssolve (const FloatColumnVector& b, octave_idx_type& info) const; FloatColumnVector lssolve (const FloatColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const; FloatColumnVector lssolve (const FloatColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const; + octave_idx_type& rank, float& rcon) const; FloatComplexColumnVector lssolve (const FloatComplexColumnVector& b) const; FloatComplexColumnVector lssolve (const FloatComplexColumnVector& b, @@ -257,7 +260,7 @@ octave_idx_type& rank) const; FloatComplexColumnVector lssolve (const FloatComplexColumnVector& b, octave_idx_type& info, - octave_idx_type& rank, float& rcond) const; + octave_idx_type& rank, float& rcon) const; FloatMatrix expm (void) const;
--- a/src/ChangeLog Tue May 13 12:40:23 2008 +0200 +++ b/src/ChangeLog Wed May 14 23:28:41 2008 +0200 @@ -1,5 +1,7 @@ 2008-05-21 David Bateman <dbateman@free.fr> + * DLD-FUNCTIONS/rcond.cc (Frcond): Add support for single precision. + * DLD-FUNCTIONS/sqrt.m: Replace DBL_* with FLT_* for single precision types. * data.cc (static octave_value fill_matrix (const
--- a/src/DLD-FUNCTIONS/rcond.cc Tue May 13 12:40:23 2008 +0200 +++ b/src/DLD-FUNCTIONS/rcond.cc Wed May 14 23:28:41 2008 +0200 @@ -51,6 +51,23 @@ 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_single_type ()) + { + if (args(0).is_complex_type ()) + { + FloatComplexMatrix m = args(0).float_complex_matrix_value (); + MatrixType mattyp; + retval = m.rcond (mattyp); + args(0).matrix_type (mattyp); + } + else + { + FloatMatrix m = args(0).float_matrix_value (); + MatrixType mattyp; + retval = m.rcond (mattyp); + args(0).matrix_type (mattyp); + } + } else if (args(0).is_complex_type ()) { ComplexMatrix m = args(0).complex_matrix_value ();