diff -r d9643bc3792a liboctave/array/CMatrix.cc --- a/liboctave/array/CMatrix.cc Thu Jul 07 19:06:09 2016 -0700 +++ b/liboctave/array/CMatrix.cc Fri Jul 08 17:29:23 2016 +0200 @@ -263,7 +263,8 @@ F77_CHAR_ARG_LEN_DECL); } -static const Complex Complex_NaN_result (octave::numeric_limits::NaN (), octave::numeric_limits::NaN ()); +static const Complex Complex_NaN_result (octave::numeric_limits::NaN (), + octave::numeric_limits::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(0)).max (); + // Calculate (always, see bug #45577) the norm of the matrix, for later use. + double anorm = + retval.abs ().sum ().row (static_cast(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::Inf (), 0.)); + ret = ComplexMatrix (rows (), columns (), + Complex (octave::numeric_limits::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 (); diff -r d9643bc3792a liboctave/array/fCMatrix.cc --- a/liboctave/array/fCMatrix.cc Thu Jul 07 19:06:09 2016 -0700 +++ b/liboctave/array/fCMatrix.cc Fri Jul 08 17:29:23 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(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(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 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(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 diff -r d9643bc3792a test/bug-46330.tst --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/bug-46330.tst Fri Jul 08 17:29:23 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 diff -r d9643bc3792a test/module.mk --- a/test/module.mk Thu Jul 07 19:06:09 2016 -0700 +++ b/test/module.mk Fri Jul 08 17:29:23 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 \