Mercurial > jwe > octave
changeset 30559:117ebe363f56 stable
Fix handling of scalar exceptional values in inv() (bug #61689)
* inv.cc (Finv): Pass "true" rather than '1' to inverse() function to match
'bool' argument. Pass 5th argument of true to inverse() to force calculation
of condition number. Rename variable "rcond_plus_one_eq_one" to "is_singular"
for clarity. Use isnan() to also catch singular matrices with NaN as
reciprocal condition number. Don't emit warning for inverse of a scalar
which thereby matches '/' and '\' operators. Add many BIST tests for various
exceptional value inputs.
* CMatrix.cc, dMatrix.cc, fCMatrix.cc, fMatrix.cc (inverse): Expand if
conditional for MatrixType::Diagonal to check input parameter "calc_cond" and
use an if/else tree to determine the reciprocal condition number (rcond) for
a scalar.
author | Rik <rik@octave.org> |
---|---|
date | Sat, 25 Dec 2021 19:16:44 -0800 |
parents | 4d2b528464ad |
children | 72cea722ca21 4054c9dc6013 |
files | libinterp/corefcn/inv.cc liboctave/array/CMatrix.cc liboctave/array/dMatrix.cc liboctave/array/fCMatrix.cc liboctave/array/fMatrix.cc |
diffstat | 5 files changed, 221 insertions(+), 58 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/inv.cc Sat Dec 25 16:10:52 2021 +0100 +++ b/libinterp/corefcn/inv.cc Sat Dec 25 19:16:44 2021 -0800 @@ -127,8 +127,8 @@ } else if (arg.is_perm_matrix ()) { + info = 0; rcond = 1.0; - info = 0; result = arg.perm_matrix_value ().inverse (); } else if (isfloat) @@ -138,7 +138,7 @@ FloatMatrix m = arg.float_matrix_value (); MatrixType mattyp = args(0).matrix_type (); - result = m.inverse (mattyp, info, frcond, 1); + result = m.inverse (mattyp, info, frcond, true, true); args(0).matrix_type (mattyp); } else if (arg.iscomplex ()) @@ -146,7 +146,7 @@ FloatComplexMatrix m = arg.float_complex_matrix_value (); MatrixType mattyp = args(0).matrix_type (); - result = m.inverse (mattyp, info, frcond, 1); + result = m.inverse (mattyp, info, frcond, true, true); args(0).matrix_type (mattyp); } } @@ -159,7 +159,7 @@ SparseMatrix m = arg.sparse_matrix_value (); MatrixType mattyp = args(0).matrix_type (); - result = m.inverse (mattyp, info, rcond, 1); + result = m.inverse (mattyp, info, rcond, true, true); args(0).matrix_type (mattyp); } else @@ -167,7 +167,7 @@ Matrix m = arg.matrix_value (); MatrixType mattyp = args(0).matrix_type (); - result = m.inverse (mattyp, info, rcond, 1); + result = m.inverse (mattyp, info, rcond, true, true); args(0).matrix_type (mattyp); } } @@ -178,7 +178,7 @@ SparseComplexMatrix m = arg.sparse_complex_matrix_value (); MatrixType mattyp = args(0).matrix_type (); - result = m.inverse (mattyp, info, rcond, 1); + result = m.inverse (mattyp, info, rcond, true, true); args(0).matrix_type (mattyp); } else @@ -186,7 +186,7 @@ ComplexMatrix m = arg.complex_matrix_value (); MatrixType mattyp = args(0).matrix_type (); - result = m.inverse (mattyp, info, rcond, 1); + result = m.inverse (mattyp, info, rcond, true, true); args(0).matrix_type (mattyp); } } @@ -200,58 +200,164 @@ if (nargout > 1) retval(1) = (isfloat ? octave_value (frcond) : octave_value (rcond)); - bool rcond_plus_one_eq_one = false; + if (nargout < 2) + { + bool is_singular; - if (isfloat) - { - volatile float xrcond = frcond; - rcond_plus_one_eq_one = xrcond + 1.0f == 1.0f; + if (isfloat) + is_singular = ((frcond + 1.0f == 1.0f) || octave::math::isnan (frcond)) + && ! arg.is_scalar_type (); + else + is_singular = ((rcond + 1.0 == 1.0) || octave::math::isnan (rcond)) + && ! arg.is_scalar_type (); + + if (info == -1 || is_singular) + warn_singular_matrix (isfloat ? frcond : rcond); } - else - { - volatile double xrcond = rcond; - rcond_plus_one_eq_one = xrcond + 1.0 == 1.0; - } - - if (nargout < 2 && (info == -1 || rcond_plus_one_eq_one)) - warn_singular_matrix (isfloat ? frcond : rcond); return retval; } /* -%!assert (inv ([1, 2; 3, 4]), [-2, 1; 1.5, -0.5], sqrt (eps)) +## Basic test for double/single matrices +%!assert (inv ([1, 2; 3, 4]), [-2, 1; 1.5, -0.5], 5*eps) +%!test +%! [xinv, rcond] = inv ([1,2;3,4]); +%! assert (xinv, [-2, 1; 1.5, -0.5], 5*eps); +%! assert (isa (rcond, "double")); + %!assert (inv (single ([1, 2; 3, 4])), single ([-2, 1; 1.5, -0.5]), -%! sqrt (eps ("single"))) +%! 5*eps ("single")) +%!test +%! [xinv, rcond] = inv (single ([1,2;3,4])); +%! assert (xinv, single ([-2, 1; 1.5, -0.5]), 5*eps ("single")); +%! assert (isa (rcond, "single")); + +## Normal scalar cases +%!assert (inv (2), 0.5) +%!test +%! [xinv, rcond] = inv (2); +%! assert (xinv, 0.5); +%! assert (rcond, 1); +%!assert (inv (single (2)), single (0.5)) +%!test +%! [xinv, rcond] = inv (single (2)); +%! assert (xinv, single (0.5)); +%! assert (rcond, single (1)); +%!assert (inv (complex (1, -1)), 0.5+0.5i) +%!test +%! [xinv, rcond] = inv (complex (1, -1)); +%! assert (xinv, 0.5+0.5i); +%! assert (rcond, 1); +%!assert (inv (complex (single (1), -1)), single (0.5+0.5i)) +%!test +%! [xinv, rcond] = inv (complex (single (1), -1)); +%! assert (xinv, single (0.5+0.5i)); +%! assert (rcond, single (1)); ## Test special inputs +## Empty matrix %!assert (inv (zeros (2,0)), []) -%!warning <matrix singular> assert (inv (0), Inf) + +## Scalar "0" +%!assert (inv (0), Inf) +%!test +%! [xinv, rcond] = inv (0); +%! assert (xinv, Inf); +%! assert (rcond, 0); +%!assert (inv (single (0)), single (Inf)) +%!test +%! [xinv, rcond] = inv (single (0)); +%! assert (xinv, single (Inf)); +%! assert (rcond, single (0)); +%!assert (inv (complex (0, 0)), Inf) +%!test +%! [xinv, rcond] = inv (complex (0, 0)); +%! assert (xinv, Inf); +%! assert (rcond, 0); +%!assert (inv (complex (single (0), 0)), single (Inf)) +%!test +%! [xinv, rcond] = inv (complex (single (0), 0)); +%! assert (xinv, single (Inf)); +%! assert (rcond, single (0)); ## NOTE: Matlab returns +Inf for -0 input, but it returns -Inf for 1/-0. -## These should be the same and in Octave they are. -%!warning <matrix singular> assert (inv (-0), -Inf) -%!warning <matrix singular> assert (inv (single (0)), single (Inf)) -%!warning <matrix singular> assert (inv (complex (0, 0)), Inf) -%!warning <matrix singular> assert (inv (single (complex (0,1)) - i), -%! single (Inf)) +## These should be the same, and in Octave they are. +%!assert (inv (-0), -Inf) +%!test +%! [xinv, rcond] = inv (-0); +%! assert (xinv, -Inf); +%! assert (rcond, 0); + +## Scalar "Inf" +%!assert (inv (Inf), 0) +%!test +%! [xinv, rcond] = inv (Inf); +%! assert (xinv, 0); +%! assert (rcond, 0); +%!assert (inv (single (Inf)), single (0)) +%!test +%! [xinv, rcond] = inv (single (Inf)); +%! assert (xinv, single (0)); +%! assert (rcond, single (0)); +%!assert (inv (complex (1, Inf)), 0) +%!test +%! [xinv, rcond] = inv (complex (1, Inf)); +%! assert (xinv, 0); +%! assert (rcond, 0); +%!assert (inv (complex (single (1), Inf)), single (0)) +%!test +%! [xinv, rcond] = inv (complex (single (1), Inf)); +%! assert (xinv, single (0)); +%! assert (rcond, single (0)); + +## Scalar "NaN" +%!assert (inv (NaN), NaN) +%!test +%! [xinv, rcond] = inv (NaN); +%! assert (xinv, NaN); +%! assert (rcond, NaN); +%!assert (inv (single (NaN)), single (NaN)) +%!test +%! [xinv, rcond] = inv (single (NaN)); +%! assert (xinv, single (NaN)); +%! assert (rcond, single (NaN)); +%!assert (inv (complex (1, NaN)), complex (NaN, NaN)) +%!test +%! [xinv, rcond] = inv (complex (1, NaN)); +%! assert (xinv, complex (NaN, NaN)); +%! assert (rcond, NaN); +%!assert (inv (complex (single (1), NaN)), complex (single (NaN), NaN)) +%!test +%! [xinv, rcond] = inv (complex (single (1), NaN)); +%! assert (xinv, complex (single (NaN), NaN)); +%! assert (rcond, single (NaN)); + +## Matrix special values +## Matrix of all zeroes %!warning <matrix singular> assert (inv (zeros (2,2)), Inf (2,2)) -%!warning <matrix singular> assert (inv (Inf), 0) -%!warning <matrix singular> assert (inv (-Inf), -0) -%!warning <matrix singular> assert (inv (single (Inf)), single (0)) -%!warning <matrix singular> assert (inv (complex (1, Inf)), 0) -%!warning <matrix singular> assert (inv (single (complex (1,Inf))), single (0)) -%!assert (inv (Inf (2,2)), NaN (2,2)) - +%!test +%! [xinv, rcond] = inv (zeros (2,2)); +%! assert (xinv, Inf (2,2)); +%! assert (rcond, 0); +## Matrix of all Inf +%!warning <rcond = > assert (inv (Inf (2,2)), NaN (2,2)) %!test -%! [xinv, rcond] = inv (single ([1,2;3,4])); -%! assert (isa (xinv, "single")); -%! assert (isa (rcond, "single")); - +%! [xinv, rcond] = inv (Inf (2,2)); +%! assert (xinv, NaN (2,2)); +%! assert (rcond, NaN); +## Matrix of all NaN +%!warning <rcond = > assert (inv (NaN (2,2)), NaN (2,2)) %!test -%! [xinv, rcond] = inv ([1,2;3,4]); -%! assert (isa (xinv, "double")); -%! assert (isa (rcond, "double")); +%! [xinv, rcond] = inv (NaN (2,2)); +%! assert (xinv, NaN (2,2)); +%! assert (rcond, NaN); +## Special diagonal matrices +%!test +%! fail ("A = inv (diag ([1, 0, 1]))", "warning", "matrix singular"); +%! assert (A, diag ([Inf, Inf, Inf])); + +## Special sparse matrices %!testif HAVE_UMFPACK <*56232> %! fail ("A = inv (sparse ([1, 2;0 ,0]))", "warning", "matrix singular"); %! assert (A, sparse ([Inf, Inf; 0, 0])); @@ -260,13 +366,6 @@ %! fail ("A = inv (sparse ([1i, 2;0 ,0]))", "warning", "matrix singular"); %! assert (A, sparse ([Inf, Inf; 0, 0])); -%!test -%! fail ("A = inv (diag ([1, 0, 1]))", "warning", "matrix singular"); -%! assert (A, diag ([Inf, Inf, Inf])); - -%!error <inverse of the null matrix not defined> inv (diag ([0, 0])) -%!error <inverse of the null matrix not defined> inv (diag (complex ([0, 0]))) - %!testif HAVE_UMFPACK <*56232> %! fail ("A = inv (sparse ([1, 0, 0; 0, 0, 0; 0, 0, 1]))", "warning", "matrix singular"); %! assert (A, sparse ([Inf, 0, 0; 0, 0, 0; 0, 0, Inf])); @@ -275,6 +374,8 @@ %!error inv ([1, 2; 3, 4], 2) %!error <must be a square matrix> inv ([1, 2; 3, 4; 5, 6]) %!error <inverse of the null matrix not defined> inv (sparse (2, 2, 0)) +%!error <inverse of the null matrix not defined> inv (diag ([0, 0])) +%!error <inverse of the null matrix not defined> inv (diag (complex ([0, 0]))) */ DEFALIAS (inverse, inv);
--- a/liboctave/array/CMatrix.cc Sat Dec 25 16:10:52 2021 +0100 +++ b/liboctave/array/CMatrix.cc Sat Dec 25 19:16:44 2021 -0800 @@ -935,13 +935,29 @@ if (typ == MatrixType::Unknown) typ = mattype.type (*this); - if (typ == MatrixType::Diagonal) // a scalar is also classified as Diagonal. + if (typ == MatrixType::Diagonal) // a scalar is classified as Diagonal. { - if (std::real (this->elem (0)) == 0 && std::imag (this->elem (0)) == 0) + Complex scalar = this->elem (0); + double real = std::real (scalar); + double imag = std::imag (scalar); + + if (real == 0 && imag == 0) ret = ComplexMatrix (1, 1, Complex (octave::numeric_limits<double>::Inf (), 0.0)); else ret = Complex (1, 0) / (*this); + + if (calc_cond) + { + if (octave::math::isfinite (real) && octave::math::isfinite (imag) + && (real != 0 || imag != 0)) + rcon = 1.0; + else if (octave::math::isinf (real) || octave::math::isinf (imag) + || (real == 0 && imag == 0)) + rcon = 0.0; + else + rcon = octave::numeric_limits<double>::NaN (); + } } else if (typ == MatrixType::Upper || typ == MatrixType::Lower) ret = tinverse (mattype, info, rcon, force, calc_cond);
--- a/liboctave/array/dMatrix.cc Sat Dec 25 16:10:52 2021 +0100 +++ b/liboctave/array/dMatrix.cc Sat Dec 25 19:16:44 2021 -0800 @@ -640,8 +640,20 @@ if (typ == MatrixType::Unknown) typ = mattype.type (*this); - if (typ == MatrixType::Diagonal) // a scalar is also classified as Diagonal. - ret = 1 / (*this); + if (typ == MatrixType::Diagonal) // a scalar is classified as Diagonal. + { + ret = 1 / (*this); + if (calc_cond) + { + double scalar = this->elem (0); + if (octave::math::isfinite (scalar) && scalar != 0) + rcon = 1.0; + else if (octave::math::isinf (scalar) || scalar == 0) + rcon = 0.0; + else + rcon = octave::numeric_limits<double>::NaN (); + } + } else if (typ == MatrixType::Upper || typ == MatrixType::Lower) ret = tinverse (mattype, info, rcon, force, calc_cond); else
--- a/liboctave/array/fCMatrix.cc Sat Dec 25 16:10:52 2021 +0100 +++ b/liboctave/array/fCMatrix.cc Sat Dec 25 19:16:44 2021 -0800 @@ -938,8 +938,30 @@ if (typ == MatrixType::Unknown) typ = mattype.type (*this); - if (typ == MatrixType::Diagonal) // a scalar is also classified as Diagonal. - ret = FloatComplex (1, 0) / (*this); + if (typ == MatrixType::Diagonal) // a scalar is classified as Diagonal. + { + FloatComplex scalar = this->elem (0); + float real = std::real (scalar); + float imag = std::imag (scalar); + + if (real == 0 && imag == 0) + ret = FloatComplexMatrix (1, 1, + FloatComplex (octave::numeric_limits<float>::Inf (), 0.0)); + else + ret = FloatComplex (1, 0) / (*this); + + if (calc_cond) + { + if (octave::math::isfinite (real) && octave::math::isfinite (imag) + && (real != 0 || imag != 0)) + rcon = 1.0f; + else if (octave::math::isinf (real) || octave::math::isinf (imag) + || (real == 0 && imag == 0)) + rcon = 0.0f; + else + rcon = octave::numeric_limits<float>::NaN (); + } + } else if (typ == MatrixType::Upper || typ == MatrixType::Lower) ret = tinverse (mattype, info, rcon, force, calc_cond); else
--- a/liboctave/array/fMatrix.cc Sat Dec 25 16:10:52 2021 +0100 +++ b/liboctave/array/fMatrix.cc Sat Dec 25 19:16:44 2021 -0800 @@ -646,8 +646,20 @@ if (typ == MatrixType::Unknown) typ = mattype.type (*this); - if (typ == MatrixType::Diagonal) // a scalar is also classified as Diagonal. - ret = 1 / (*this); + if (typ == MatrixType::Diagonal) // a scalar is classified as Diagonal. + { + ret = 1 / (*this); + if (calc_cond) + { + float scalar = this->elem (0); + if (octave::math::isfinite (scalar) && scalar != 0) + rcon = 1.0f; + else if (octave::math::isinf (scalar) || scalar == 0) + rcon = 0.0f; + else + rcon = octave::numeric_limits<float>::NaN (); + } + } else if (typ == MatrixType::Upper || typ == MatrixType::Lower) ret = tinverse (mattype, info, rcon, force, calc_cond); else