Mercurial > octave

changeset 30554:117ebe363f56 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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