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;