# HG changeset patch # User Rik # Date 1637542429 28800 # Node ID 98400baa509fc5c0a25a05bebb25621b342eb66b # Parent 90ca4cd8e7e5e833ff4c1e9d66db221fc7610d31 Return Inf for inv(0) for compatibility with division operator and Matlab (bug #52285). * CMatrix.cc, dMatrix.cc, fCMatrix.cc, fMatrix.cc (inverse): Check early for MatrixType::Diagonal which implies a scalar. In scalar case, just calculate 1 / *this. * inv.cc (Finverse): Add BIST tests for new behavior. diff -r 90ca4cd8e7e5 -r 98400baa509f libinterp/corefcn/inv.cc --- a/libinterp/corefcn/inv.cc Sun Nov 21 13:39:19 2021 -0800 +++ b/libinterp/corefcn/inv.cc Sun Nov 21 16:53:49 2021 -0800 @@ -221,15 +221,26 @@ /* %!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"))) +%!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 assert (inv (0), Inf) +## 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 assert (inv (-0), -Inf) +%!warning assert (inv (single (0)), single (Inf)) +%!warning assert (inv (complex (0, 0)), Inf) +%!warning assert (inv (single (complex (0,1)) - i), +%! single (Inf)) +%!warning assert (inv (zeros (2,2)), Inf (2,2)) %!warning assert (inv (Inf), 0) %!warning assert (inv (-Inf), -0) %!warning assert (inv (single (Inf)), single (0)) %!warning assert (inv (complex (1, Inf)), 0) %!warning assert (inv (single (complex (1,Inf))), single (0)) +%!assert (inv (Inf (2,2)), NaN (2,2)) %!test %! [xinv, rcond] = inv (single ([1,2;3,4])); diff -r 90ca4cd8e7e5 -r 98400baa509f liboctave/array/CMatrix.cc --- a/liboctave/array/CMatrix.cc Sun Nov 21 13:39:19 2021 -0800 +++ b/liboctave/array/CMatrix.cc Sun Nov 21 16:53:49 2021 -0800 @@ -935,7 +935,15 @@ if (typ == MatrixType::Unknown) typ = mattype.type (*this); - if (typ == MatrixType::Upper || typ == MatrixType::Lower) + if (typ == MatrixType::Diagonal) // a scalar is also classified as Diagonal. + { + if (std::real (this->elem (0)) == 0 && std::imag (this->elem (0)) == 0) + ret = ComplexMatrix (1,1, + Complex (octave::numeric_limits::Inf (), 0.0)); + else + ret = Complex (1,0) / (*this); + } + else if (typ == MatrixType::Upper || typ == MatrixType::Lower) ret = tinverse (mattype, info, rcon, force, calc_cond); else { @@ -959,11 +967,8 @@ if ((calc_cond || mattype.ishermitian ()) && rcon == 0.0) { - if (numel () == 1) - ret = ComplexMatrix (1, 1, 0.0); - else - ret = ComplexMatrix (rows (), columns (), - Complex (octave::numeric_limits::Inf (), 0.0)); + ret = ComplexMatrix (rows (), columns (), + Complex (octave::numeric_limits::Inf (), 0.0)); } } diff -r 90ca4cd8e7e5 -r 98400baa509f liboctave/array/dMatrix.cc --- a/liboctave/array/dMatrix.cc Sun Nov 21 13:39:19 2021 -0800 +++ b/liboctave/array/dMatrix.cc Sun Nov 21 16:53:49 2021 -0800 @@ -640,7 +640,9 @@ if (typ == MatrixType::Unknown) typ = mattype.type (*this); - if (typ == MatrixType::Upper || typ == MatrixType::Lower) + if (typ == MatrixType::Diagonal) // a scalar is also classified as Diagonal. + ret = 1 / (*this); + else if (typ == MatrixType::Upper || typ == MatrixType::Lower) ret = tinverse (mattype, info, rcon, force, calc_cond); else { @@ -662,8 +664,7 @@ if (! mattype.ishermitian ()) ret = finverse (mattype, info, rcon, force, calc_cond); - if ((calc_cond || mattype.ishermitian ()) && rcon == 0.0 - && (numel () != 1)) + if ((calc_cond || mattype.ishermitian ()) && rcon == 0.0) ret = Matrix (rows (), columns (), octave::numeric_limits::Inf ()); } diff -r 90ca4cd8e7e5 -r 98400baa509f liboctave/array/fCMatrix.cc --- a/liboctave/array/fCMatrix.cc Sun Nov 21 13:39:19 2021 -0800 +++ b/liboctave/array/fCMatrix.cc Sun Nov 21 16:53:49 2021 -0800 @@ -938,7 +938,9 @@ if (typ == MatrixType::Unknown) typ = mattype.type (*this); - if (typ == MatrixType::Upper || typ == MatrixType::Lower) + if (typ == MatrixType::Diagonal) // a scalar is also classified as Diagonal. + ret = FloatComplex (1,0) / (*this); + else if (typ == MatrixType::Upper || typ == MatrixType::Lower) ret = tinverse (mattype, info, rcon, force, calc_cond); else { @@ -962,11 +964,8 @@ if ((calc_cond || mattype.ishermitian ()) && rcon == 0.0) { - if (numel () == 1) - ret = FloatComplexMatrix (1, 1, 0.0); - else - ret = FloatComplexMatrix (rows (), columns (), - FloatComplex (octave::numeric_limits::Inf (), 0.0)); + ret = FloatComplexMatrix (rows (), columns (), + FloatComplex (octave::numeric_limits::Inf (), 0.0)); } } diff -r 90ca4cd8e7e5 -r 98400baa509f liboctave/array/fMatrix.cc --- a/liboctave/array/fMatrix.cc Sun Nov 21 13:39:19 2021 -0800 +++ b/liboctave/array/fMatrix.cc Sun Nov 21 16:53:49 2021 -0800 @@ -646,7 +646,9 @@ if (typ == MatrixType::Unknown) typ = mattype.type (*this); - if (typ == MatrixType::Upper || typ == MatrixType::Lower) + if (typ == MatrixType::Diagonal) // a scalar is also classified as Diagonal. + ret = 1 / (*this); + else if (typ == MatrixType::Upper || typ == MatrixType::Lower) ret = tinverse (mattype, info, rcon, force, calc_cond); else { @@ -668,8 +670,7 @@ if (! mattype.ishermitian ()) ret = finverse (mattype, info, rcon, force, calc_cond); - if ((calc_cond || mattype.ishermitian ()) && rcon == 0.0 - && (numel () != 1)) + if ((calc_cond || mattype.ishermitian ()) && rcon == 0.0) ret = FloatMatrix (rows (), columns (), octave::numeric_limits::Inf ()); }