Mercurial > octave

changeset 22085:bf5fbf347aaf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Avoid segfault with complex matrices using LAPACK (bug #46330). * liboctave/array/CMatrix.cc: code refactoring. * liboctave/array/fCMatrix.cc(FloatComplexMatrix::determinant, FloatComplexMatrix::rcond): compare Matrix norm with Nan. * liboctave/array/fCMatrix.cc(FloatComplexMatrix::finverse, FloatComplexMatrix::fsolve): compare Matrix norm with Inf and Nan. * test/module.mk: Add new test file. * test/bug-46330.tst: Add regression test case from bug report. For Windows this test will be omitted (see bug #39000, Can't override BLAS XERBLA handler).
author Kai T. Ohlhus <k.ohlhus@gmail.com>
date Mon, 11 Jul 2016 09:27:10 +0200
parents 79ee6df71b51
children 67a44207da71
files liboctave/array/CMatrix.cc liboctave/array/fCMatrix.cc test/bug-46330.tst test/module.mk
diffstat 4 files changed, 78 insertions(+), 86 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/array/CMatrix.cc	Sun Jul 10 11:02:44 2016 -0700
+++ b/liboctave/array/CMatrix.cc	Mon Jul 11 09:27:10 2016 +0200
@@ -263,7 +263,8 @@
                                F77_CHAR_ARG_LEN_DECL);
 }
 
-static const Complex Complex_NaN_result (octave::numeric_limits<double>::NaN (), octave::numeric_limits<double>::NaN ());
+static const Complex Complex_NaN_result (octave::numeric_limits<double>::NaN (),
+                                         octave::numeric_limits<double>::NaN ());
 
 // Complex Matrix class
 
@@ -1043,10 +1044,9 @@
 
   info = 0;
 
-  // Calculate the norm of the matrix, for later use.
-  double anorm;
-  //if (calc_cond)   // Must always calculate anorm for bug #45577
-  anorm = retval.abs ().sum ().row (static_cast<octave_idx_type>(0)).max ();
+  // Calculate (always, see bug #45577) the norm of the matrix, for later use.
+  double anorm =
+    retval.abs ().sum ().row (static_cast<octave_idx_type>(0)).max ();
 
   // Work around bug #45577, LAPACK crashes Octave if norm is NaN
   // and bug #46330, segfault with matrices containing Inf & NaN
@@ -1127,7 +1127,8 @@
         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.));
+        ret = ComplexMatrix (rows (), columns (),
+                             Complex (octave::numeric_limits<double>::Inf (), 0.));
     }
 
   return ret;
@@ -1611,10 +1612,8 @@
 
       info = 0;
 
-      // Calculate the norm of the matrix, for later use.
-      double anorm = 0;
-      //if (calc_cond)   // Must always calculate anorm for bug #45577
-      anorm = xnorm (*this, 1);
+      // Calculate (always, see bug #45577) the norm of the matrix, for later use.
+      double anorm = xnorm (*this, 1);
 
       // Work around bug #45577, LAPACK crashes Octave if norm is NaN
       if (octave::math::isnan (anorm))
@@ -2203,6 +2202,7 @@
                 mattype.mark_as_rectangular ();
             }
         }
+
       if (octave::math::isinf (anorm))
         {
           retval = ComplexMatrix (b.rows (), b.cols (), Complex (0, 0));
@@ -3733,21 +3733,7 @@
 ComplexMatrix
 min (const ComplexMatrix& m, const Complex& c)
 {
-  octave_idx_type nr = m.rows ();
-  octave_idx_type nc = m.columns ();
-
-  EMPTY_RETURN_CHECK (ComplexMatrix);
-
-  ComplexMatrix result (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    for (octave_idx_type i = 0; i < nr; i++)
-      {
-        octave_quit ();
-        result(i, j) = octave::math::min (m(i, j), c);
-      }
-
-  return result;
+  return min (c, m);
 }
 
 ComplexMatrix
@@ -3818,21 +3804,7 @@
 ComplexMatrix
 max (const ComplexMatrix& m, const Complex& c)
 {
-  octave_idx_type nr = m.rows ();
-  octave_idx_type nc = m.columns ();
-
-  EMPTY_RETURN_CHECK (ComplexMatrix);
-
-  ComplexMatrix result (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    for (octave_idx_type i = 0; i < nr; i++)
-      {
-        octave_quit ();
-        result(i, j) = octave::math::max (m(i, j), c);
-      }
-
-  return result;
+  return max (c, m);
 }
 
 ComplexMatrix
@@ -3862,6 +3834,7 @@
             }
         }
 
+      // FIXME: is it so much faster?
       if (columns_are_real_only)
         {
           for (octave_idx_type i = 0; i < nr; i++)
@@ -3886,7 +3859,6 @@
 ComplexMatrix linspace (const ComplexColumnVector& x1,
                         const ComplexColumnVector& x2,
                         octave_idx_type n)
-
 {
   octave_idx_type m = x1.numel ();
 
--- a/liboctave/array/fCMatrix.cc	Sun Jul 10 11:02:44 2016 -0700
+++ b/liboctave/array/fCMatrix.cc	Mon Jul 11 09:27:10 2016 +0200
@@ -1049,13 +1049,16 @@
 
   info = 0;
 
-  // Calculate the norm of the matrix, for later use.
-  float anorm;
-  if (calc_cond)
-    anorm = retval.abs ().sum ().row (static_cast<octave_idx_type>(0))
-            .max ();
-
-  F77_XFCN (cgetrf, CGETRF, (nc, nc, tmp_data, nr, pipvt, info));
+  // Calculate (always, see bug #45577) the norm of the matrix, for later use.
+  float anorm =
+    retval.abs ().sum ().row (static_cast<octave_idx_type>(0)).max ();
+
+  // Work around bug #45577, LAPACK crashes Octave if norm is NaN
+  // and bug #46330, segfault with matrices containing Inf & NaN
+  if (octave::math::isnan (anorm) || octave::math::isinf (anorm))
+    info = -1;
+  else
+    F77_XFCN (cgetrf, CGETRF, (nc, nc, tmp_data, nr, pipvt, info));
 
   // Throw-away extra info LAPACK gives so as to not change output.
   rcon = 0.0;
@@ -1077,7 +1080,7 @@
         info = -1;
     }
 
-  if (info == -1 && ! force)
+  if ((info == -1 && ! force) || octave::math::isinf (anorm))
     retval = *this;  // Restore contents.
   else
     {
@@ -1611,11 +1614,14 @@
 
       info = 0;
 
-      // Calculate the norm of the matrix, for later use.
-      float anorm = 0;
-      if (calc_cond) anorm = xnorm (*this, 1);
-
-      F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info));
+      // Calculate (always, see bug #45577) the norm of the matrix, for later use.
+      float anorm = xnorm (*this, 1);
+
+      // Work around bug #45577, LAPACK crashes Octave if norm is NaN
+      if (octave::math::isnan (anorm))
+        info = -1;
+      else
+        F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info));
 
       // Throw-away extra info LAPACK gives so as to not change output.
       rcon = 0.0;
@@ -1805,7 +1811,11 @@
               Array<float> rz (dim_vector (2 * nc, 1));
               float *prz = rz.fortran_vec ();
 
-              F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info));
+              // Work around bug #45577, LAPACK crashes Octave if norm is NaN
+              if (octave::math::isnan (anorm))
+                info = -1;
+              else
+                F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info));
 
               if (info != 0)
                 {
@@ -2137,7 +2147,12 @@
             anorm = atmp.abs ().sum ().row (static_cast<octave_idx_type>(0))
                     .max ();
 
-          F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info));
+          // Work around bug #45577, LAPACK crashes Octave if norm is NaN
+          // and bug #46330, segfault with matrices containing Inf & NaN
+          if (octave::math::isnan (anorm) || octave::math::isinf (anorm))
+            info = -2;
+          else
+            F77_XFCN (cgetrf, CGETRF, (nr, nr, tmp_data, nr, pipvt, info));
 
           // Throw-away extra info LAPACK gives so as to not change output.
           rcon = 0.0;
@@ -2197,6 +2212,13 @@
                 mattype.mark_as_rectangular ();
             }
         }
+
+      if (octave::math::isinf (anorm))
+        {
+          retval = FloatComplexMatrix (b.rows (), b.cols (),
+                                       FloatComplex (0, 0));
+          mattype.mark_as_full ();
+        }
     }
 
   return retval;
@@ -3742,21 +3764,7 @@
 FloatComplexMatrix
 min (const FloatComplexMatrix& m, const FloatComplex& c)
 {
-  octave_idx_type nr = m.rows ();
-  octave_idx_type nc = m.columns ();
-
-  EMPTY_RETURN_CHECK (FloatComplexMatrix);
-
-  FloatComplexMatrix result (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    for (octave_idx_type i = 0; i < nr; i++)
-      {
-        octave_quit ();
-        result(i, j) = octave::math::min (m(i, j), c);
-      }
-
-  return result;
+  return min (c, m);
 }
 
 FloatComplexMatrix
@@ -3827,21 +3835,7 @@
 FloatComplexMatrix
 max (const FloatComplexMatrix& m, const FloatComplex& c)
 {
-  octave_idx_type nr = m.rows ();
-  octave_idx_type nc = m.columns ();
-
-  EMPTY_RETURN_CHECK (FloatComplexMatrix);
-
-  FloatComplexMatrix result (nr, nc);
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    for (octave_idx_type i = 0; i < nr; i++)
-      {
-        octave_quit ();
-        result(i, j) = octave::math::max (m(i, j), c);
-      }
-
-  return result;
+  return max (c, m);
 }
 
 FloatComplexMatrix
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/bug-46330.tst	Mon Jul 11 09:27:10 2016 +0200
@@ -0,0 +1,25 @@
+## bug #46330: segfault with matrices containing Inf & NaN
+
+%!warning
+%! a = [-0.46E-01,            0.10681415316, 0.0,   -0.17121680433;
+%!      -0.1675901504661613, -0.515,         1.0,    0.6420630320636088E-02;
+%!       0.1543104215347786, -0.547945,     -0.906, -0.1521689385990753E-02;
+%!       0.0,                 0.0,           1.0,    0.0];
+%!
+%! b = [0.1602300107479095,      0.2111848453E-02;
+%!      0.8196877780963616E-02, -0.3025E-01;
+%!      0.9173594317692437E-01, -0.75283075;
+%!      0.0,                     0.0];
+%!
+%! c = [1.0, 0.0, 0.0, 0.0;
+%!      0.0, 0.0, 0.0, 1.0];
+%!
+%! ## This statement caused an error in LAPACK and eventually caused
+%! ## a segmentation fault.
+%! ## Triggers "warning: matrix singular to machine precision"
+%! ## FIXME: LAPACK errors become fatal crashes on Windows, don't test this
+%! if (ispc ())
+%!   warning ("unable to test for bug #46330 on Windows");
+%! else
+%!   assert (c / (i * Inf * eye (4) - a) * b, zeros (2, 2))
+%! endif
--- a/test/module.mk	Sun Jul 10 11:02:44 2016 -0700
+++ b/test/module.mk	Mon Jul 11 09:27:10 2016 +0200
@@ -10,6 +10,7 @@
   test/bug-31371.tst \
   test/bug-38565.tst \
   test/bug-38576.tst \
+  test/bug-46330.tst \
   test/colormaps.tst \
   test/command.tst \
   test/complex.tst \