Mercurial > octave
changeset 22849:1d242ae72240
Return 0 for special case of scalar Inf input to inverse (bug #49690).
* inv.cc (Finv): Add BIST tests for new behavior.
* CMatrix.cc, fCMatrix.cc, dMatrix.cc, fMatrix.cc (::inverse):
If function has had to calculate the reciprocal condition number, and
the rcond value is zero, then return 0 if the input was a scalar, or
a matrix of all Inf values.
author | Rik <rik@octave.org> |
---|---|
date | Tue, 29 Nov 2016 16:46:35 -0800 |
parents | a3bd04422a87 |
children | 76fcce30dd32 |
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, 42 insertions(+), 19 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/inv.cc Tue Nov 29 16:07:05 2016 -0500 +++ b/libinterp/corefcn/inv.cc Tue Nov 29 16:46:35 2016 -0800 @@ -207,19 +207,28 @@ %!assert (inv ([1, 2; 3, 4]), [-2, 1; 1.5, -0.5], sqrt (eps)) %!assert (inv (single ([1, 2; 3, 4])), single ([-2, 1; 1.5, -0.5]), sqrt (eps ("single"))) +## Test special inputs +%!assert (inv (zeros (2,0)), []) +%!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)) + +%!test +%! [xinv, rcond] = inv (single ([1,2;3,4])); +%! assert (isa (xinv, "single")); +%! assert (isa (rcond, "single")); + +%!test +%! [xinv, rcond] = inv ([1,2;3,4]); +%! assert (isa (xinv, "double")); +%! assert (isa (rcond, "double")); + %!error inv () %!error inv ([1, 2; 3, 4], 2) %!error <must be a square matrix> inv ([1, 2; 3, 4; 5, 6]) -%!test -%! [xinv, rcond] = inv (single ([1,2;3,4])); -%! assert (isa (xinv, 'single')); -%! assert (isa (rcond, 'single')); - -%!test -%! [xinv, rcond] = inv ([1,2;3,4]); -%! assert (isa (xinv, 'double')); -%! assert (isa (rcond, 'double')); */ // FIXME: this should really be done with an alias, but
--- a/liboctave/array/CMatrix.cc Tue Nov 29 16:07:05 2016 -0500 +++ b/liboctave/array/CMatrix.cc Tue Nov 29 16:46:35 2016 -0800 @@ -924,9 +924,14 @@ if (! mattype.is_hermitian ()) ret = finverse (mattype, info, rcon, force, calc_cond); - if ((mattype.is_hermitian () || calc_cond) && rcon == 0.) - ret = ComplexMatrix (rows (), columns (), - Complex (octave::numeric_limits<double>::Inf (), 0.)); + if ((calc_cond || mattype.is_hermitian ()) && rcon == 0.) + { + if (numel () == 1) + ret = ComplexMatrix (1, 1, 0.); + else + ret = ComplexMatrix (rows (), columns (), + Complex (octave::numeric_limits<double>::Inf (), 0.)); + } } return ret;
--- a/liboctave/array/dMatrix.cc Tue Nov 29 16:07:05 2016 -0500 +++ b/liboctave/array/dMatrix.cc Tue Nov 29 16:46:35 2016 -0800 @@ -632,8 +632,10 @@ if (! mattype.is_hermitian ()) ret = finverse (mattype, info, rcon, force, calc_cond); - if ((mattype.is_hermitian () || calc_cond) && rcon == 0.) - ret = Matrix (rows (), columns (), octave::numeric_limits<double>::Inf ()); + if ((calc_cond || mattype.is_hermitian ()) && rcon == 0. + && ! (numel () == 1)) + ret = Matrix (rows (), columns (), + octave::numeric_limits<double>::Inf ()); } return ret;
--- a/liboctave/array/fCMatrix.cc Tue Nov 29 16:07:05 2016 -0500 +++ b/liboctave/array/fCMatrix.cc Tue Nov 29 16:46:35 2016 -0800 @@ -926,9 +926,14 @@ if (! mattype.is_hermitian ()) ret = finverse (mattype, info, rcon, force, calc_cond); - if ((mattype.is_hermitian () || calc_cond) && rcon == 0.) - ret = FloatComplexMatrix (rows (), columns (), - FloatComplex (octave::numeric_limits<float>::Inf (), 0.)); + if ((calc_cond || mattype.is_hermitian ()) && rcon == 0.) + { + if (numel () == 1) + ret = FloatComplexMatrix (1, 1, 0.); + else + ret = FloatComplexMatrix (rows (), columns (), + FloatComplex (octave::numeric_limits<float>::Inf (), 0.)); + } } return ret;
--- a/liboctave/array/fMatrix.cc Tue Nov 29 16:07:05 2016 -0500 +++ b/liboctave/array/fMatrix.cc Tue Nov 29 16:46:35 2016 -0800 @@ -638,8 +638,10 @@ if (! mattype.is_hermitian ()) ret = finverse (mattype, info, rcon, force, calc_cond); - if ((mattype.is_hermitian () || calc_cond) && rcon == 0.) - ret = FloatMatrix (rows (), columns (), octave::numeric_limits<float>::Inf ()); + if ((calc_cond || mattype.is_hermitian ()) && rcon == 0. + && ! (numel () == 1)) + ret = FloatMatrix (rows (), columns (), + octave::numeric_limits<float>::Inf ()); } return ret;