Mercurial > octave

changeset 30313:98400baa509f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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.
author Rik <rik@octave.org>
date Sun, 21 Nov 2021 16:53:49 -0800
parents 90ca4cd8e7e5
children 3e09b065779d
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, 36 insertions(+), 19 deletions(-) [+]
line wrap: on
line diff
--- 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 <matrix singular> 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 <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))
+%!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 (single ([1,2;3,4]));
--- 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<double>::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<double>::Inf (), 0.0));
+          ret = ComplexMatrix (rows (), columns (),
+                               Complex (octave::numeric_limits<double>::Inf (), 0.0));
         }
     }
 
--- 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<double>::Inf ());
     }
--- 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<float>::Inf (), 0.0));
+          ret = FloatComplexMatrix (rows (), columns (),
+                                    FloatComplex (octave::numeric_limits<float>::Inf (), 0.0));
         }
     }
 
--- 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<float>::Inf ());
     }