Mercurial > hg > octave-shane
changeset 15382:197774b411ec stable
rcond: use new copy of data for full factorization if positive definite cholesky factorization fails (bug #37336)
* dMatrix.cc (Matrix::rcond): Don't reuse modified matrix data if
positive definite cholesky factorization was attempted but fails.
* CMatrix.cc (ComplexMatrix::rcond): Likewise.
* fMatrix.cc (FloatMatrix::rcond): Likewise.
* fCMatrix.cc (FloatComplexMatrix::rcond): Likewise.
* rcond.cc: New tests.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Thu, 13 Sep 2012 15:14:47 -0400 |
parents | 4663cc835c65 |
children | 2137f5638521 fd5c0159b588 |
files | liboctave/CMatrix.cc liboctave/dMatrix.cc liboctave/fCMatrix.cc liboctave/fMatrix.cc src/DLD-FUNCTIONS/rcond.cc |
diffstat | 5 files changed, 46 insertions(+), 8 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -1786,13 +1786,15 @@ 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'; + + ComplexMatrix atmp = *this; + Complex *tmp_data = atmp.fortran_vec (); + anorm = atmp.abs().sum(). row(static_cast<octave_idx_type>(0)).max(); @@ -1829,6 +1831,9 @@ { octave_idx_type info = 0; + ComplexMatrix atmp = *this; + Complex *tmp_data = atmp.fortran_vec (); + Array<octave_idx_type> ipvt (dim_vector (nr, 1)); octave_idx_type *pipvt = ipvt.fortran_vec (); @@ -2098,8 +2103,10 @@ { info = 0; char job = 'L'; + ComplexMatrix atmp = *this; Complex *tmp_data = atmp.fortran_vec (); + anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -1454,13 +1454,15 @@ 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'; + + Matrix atmp = *this; + double *tmp_data = atmp.fortran_vec (); + anorm = atmp.abs().sum(). row(static_cast<octave_idx_type>(0)).max(); @@ -1495,6 +1497,9 @@ { octave_idx_type info = 0; + Matrix atmp = *this; + double *tmp_data = atmp.fortran_vec (); + Array<octave_idx_type> ipvt (dim_vector (nr, 1)); octave_idx_type *pipvt = ipvt.fortran_vec (); @@ -1760,8 +1765,10 @@ { info = 0; char job = 'L'; + Matrix atmp = *this; double *tmp_data = atmp.fortran_vec (); + anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, @@ -1838,6 +1845,7 @@ Matrix atmp = *this; double *tmp_data = atmp.fortran_vec (); + if(anorm < 0.) anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
--- a/liboctave/fCMatrix.cc +++ b/liboctave/fCMatrix.cc @@ -1782,13 +1782,15 @@ 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'; + + FloatComplexMatrix atmp = *this; + FloatComplex *tmp_data = atmp.fortran_vec (); + anorm = atmp.abs().sum(). row(static_cast<octave_idx_type>(0)).max(); @@ -1825,6 +1827,9 @@ { octave_idx_type info = 0; + FloatComplexMatrix atmp = *this; + FloatComplex *tmp_data = atmp.fortran_vec (); + Array<octave_idx_type> ipvt (dim_vector (nr, 1)); octave_idx_type *pipvt = ipvt.fortran_vec (); @@ -2094,8 +2099,10 @@ { info = 0; char job = 'L'; + FloatComplexMatrix atmp = *this; FloatComplex *tmp_data = atmp.fortran_vec (); + anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
--- a/liboctave/fMatrix.cc +++ b/liboctave/fMatrix.cc @@ -1454,13 +1454,15 @@ 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'; + + FloatMatrix atmp = *this; + float *tmp_data = atmp.fortran_vec (); + anorm = atmp.abs().sum(). row(static_cast<octave_idx_type>(0)).max(); @@ -1495,6 +1497,9 @@ { octave_idx_type info = 0; + FloatMatrix atmp = *this; + float *tmp_data = atmp.fortran_vec (); + Array<octave_idx_type> ipvt (dim_vector (nr, 1)); octave_idx_type *pipvt = ipvt.fortran_vec (); @@ -1760,8 +1765,10 @@ { info = 0; char job = 'L'; + FloatMatrix atmp = *this; float *tmp_data = atmp.fortran_vec (); + anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max(); F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, @@ -1838,6 +1845,7 @@ FloatMatrix atmp = *this; float *tmp_data = atmp.fortran_vec (); + if(anorm < 0.) anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
--- a/src/DLD-FUNCTIONS/rcond.cc +++ b/src/DLD-FUNCTIONS/rcond.cc @@ -93,4 +93,12 @@ %!assert( rcond ([1 1; 2 1]), 1/9) %!assert( rcond (magic (4)), 0, eps) +%!shared x, sx +%! x = [-5.25, -2.25; -2.25, 1] * eps () + ones (2) / 2; +%! sx = [-5.25, -2.25; -2.25, 1] * eps ("single") + ones (2) / 2; +%!assert (rcond (x) < eps ()); +%!assert (rcond (sx) < eps ('single')); +%!assert (rcond (x*i) < eps ()); +%!assert (rcond (sx*i) < eps ('single')); + */