Mercurial > octave

changeset 24311:0643533930e7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Return NaN arrays for left division operator when operands contain NaNs (bug #52516). * CMatrix.cc, dMatrix.cc, fCMatrix.cc, fMatrix.cc (norm1): Recode to shortircuit on either NaN or Inf and return norm . * CMatrix.cc, dMatrix.cc, fCMatrix.cc, fMatrix.cc (lssolve): Check norm of matrix and return NaN matrix for NaN norm and zeros matrix for Inf norm.
author Rik <rik@octave.org>
date Sun, 26 Nov 2017 14:50:20 -0800
parents 4dce9d03b2ba
children 9d25e88d83f6
files liboctave/array/CMatrix.cc liboctave/array/dMatrix.cc liboctave/array/fCMatrix.cc liboctave/array/fMatrix.cc
diffstat 4 files changed, 50 insertions(+), 23 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/array/CMatrix.cc	Sun Nov 26 10:22:39 2017 -0800
+++ b/liboctave/array/CMatrix.cc	Sun Nov 26 14:50:20 2017 -0800
@@ -714,18 +714,19 @@
 double
 norm1 (const ComplexMatrix& a)
 {
+  double anorm = 0.0;
   ColumnVector colsum = a.abs ().sum ().row (0);
-  double anorm = -octave::numeric_limits<double>::Inf ();
 
   for (octave_idx_type i = 0; i < colsum.numel (); i++)
     {
-      if (octave::math::isnan (colsum.xelem (i)))
+      double sum = colsum.xelem (i);
+      if (octave::math::isinf (sum) || octave::math::isnan (sum))
         {
-          anorm = octave::numeric_limits<double>::NaN ();
+          anorm = sum;  // Pass Inf or NaN to output
           break;
         }
       else
-        anorm = std::max (anorm, colsum.xelem (i));
+        anorm = std::max (anorm, sum);
     }
 
   return anorm;
@@ -2574,11 +2575,17 @@
 
       anorm = norm1 (*this);
 
-      if (octave::math::isinf (anorm) || octave::math::isnan (anorm))
+      if (octave::math::isinf (anorm))
         {
           rcon = 0.0;
           octave::warn_singular_matrix ();
-          retval = Matrix (n, b_nc, 0.0);
+          retval = ComplexMatrix (n, b_nc, 0.0);
+        }
+      else if (octave::math::isnan (anorm))
+        {
+          rcon = octave::numeric_limits<double>::NaN ();
+          retval = ComplexMatrix (n, b_nc,
+                                  octave::numeric_limits<double>::NaN ());
         }
       else
         {
--- a/liboctave/array/dMatrix.cc	Sun Nov 26 10:22:39 2017 -0800
+++ b/liboctave/array/dMatrix.cc	Sun Nov 26 14:50:20 2017 -0800
@@ -428,18 +428,19 @@
 double
 norm1 (const Matrix& a)
 {
+  double anorm = 0.0;
   ColumnVector colsum = a.abs ().sum ().row (0);
-  double anorm = -octave::numeric_limits<double>::Inf ();
 
   for (octave_idx_type i = 0; i < colsum.numel (); i++)
     {
-      if (octave::math::isnan (colsum.xelem (i)))
+      double sum = colsum.xelem (i);
+      if (octave::math::isinf (sum) || octave::math::isnan (sum))
         {
-          anorm = octave::numeric_limits<double>::NaN ();
+          anorm = sum;  // Pass Inf or NaN to output
           break;
         }
       else
-        anorm = std::max (anorm, colsum.xelem (i));
+        anorm = std::max (anorm, sum);
     }
 
   return anorm;
@@ -2226,12 +2227,17 @@
 
       anorm = norm1 (*this);
 
-      if (octave::math::isinf (anorm) || octave::math::isnan (anorm))
+      if (octave::math::isinf (anorm))
         {
           rcon = 0.0;
           octave::warn_singular_matrix ();
           retval = Matrix (n, b_nc, 0.0);
         }
+      else if (octave::math::isnan (anorm))
+        {
+          rcon = octave::numeric_limits<double>::NaN ();
+          retval = Matrix (n, b_nc, octave::numeric_limits<double>::NaN ());
+        }
       else
         {
           F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval,
--- a/liboctave/array/fCMatrix.cc	Sun Nov 26 10:22:39 2017 -0800
+++ b/liboctave/array/fCMatrix.cc	Sun Nov 26 14:50:20 2017 -0800
@@ -717,18 +717,19 @@
 float
 norm1 (const FloatComplexMatrix& a)
 {
+  float anorm = 0.0;
   FloatColumnVector colsum = a.abs ().sum ().row (0);
-  float anorm = -octave::numeric_limits<float>::Inf ();
 
   for (octave_idx_type i = 0; i < colsum.numel (); i++)
     {
-      if (octave::math::isnan (colsum.xelem (i)))
+      float sum = colsum.xelem (i);
+      if (octave::math::isinf (sum) || octave::math::isnan (sum))
         {
-          anorm = octave::numeric_limits<float>::NaN ();
+          anorm = sum;  // Pass Inf or NaN to output
           break;
         }
       else
-        anorm = std::max (anorm, colsum.xelem (i));
+        anorm = std::max (anorm, sum);
     }
 
   return anorm;
@@ -2595,11 +2596,17 @@
 
       anorm = norm1 (*this);
 
-      if (octave::math::isinf (anorm) || octave::math::isnan (anorm))
+      if (octave::math::isinf (anorm))
         {
           rcon = 0.0;
           octave::warn_singular_matrix ();
-          retval = Matrix (n, b_nc, 0.0);
+          retval = FloatComplexMatrix (n, b_nc, 0.0);
+        }
+      else if (octave::math::isnan (anorm))
+        {
+          rcon = octave::numeric_limits<float>::NaN ();
+          retval = FloatComplexMatrix (n, b_nc,
+                                       octave::numeric_limits<float>::NaN ());
         }
       else
         {
--- a/liboctave/array/fMatrix.cc	Sun Nov 26 10:22:39 2017 -0800
+++ b/liboctave/array/fMatrix.cc	Sun Nov 26 14:50:20 2017 -0800
@@ -434,18 +434,19 @@
 float
 norm1 (const FloatMatrix& a)
 {
+  float anorm = 0.0;
   FloatColumnVector colsum = a.abs ().sum ().row (0);
-  float anorm = -octave::numeric_limits<float>::Inf ();
 
   for (octave_idx_type i = 0; i < colsum.numel (); i++)
     {
-      if (octave::math::isnan (colsum.xelem (i)))
+      float sum = colsum.xelem (i);
+      if (octave::math::isinf (sum) || octave::math::isnan (sum))
         {
-          anorm = octave::numeric_limits<float>::NaN ();
+          anorm = sum;  // Pass Inf or NaN to output
           break;
         }
       else
-        anorm = std::max (anorm, colsum.xelem (i));
+        anorm = std::max (anorm, sum);
     }
 
   return anorm;
@@ -2251,11 +2252,17 @@
 
       anorm = norm1 (*this);
 
-      if (octave::math::isinf (anorm) || octave::math::isnan (anorm))
+      if (octave::math::isinf (anorm))
         {
           rcon = 0.0;
           octave::warn_singular_matrix ();
-          retval = Matrix (n, b_nc, 0.0);
+          retval = FloatMatrix (n, b_nc, 0.0);
+        }
+      else if (octave::math::isnan (anorm))
+        {
+          rcon = octave::numeric_limits<float>::NaN ();
+          retval = FloatMatrix (n, b_nc,
+                                octave::numeric_limits<float>::NaN ());
         }
       else
         {