Mercurial > octave
changeset 22598:5aa8f199e328 stable
avoid invalid BLAS calls that then invoke xerbla (bug #39000)
* dMatrix.cc (Matrix::lssolve): Don't call dgelsd if norm is inf or nan.
* fMatrix.cc (FloatMatrix::lssolve): Likewise, for sgelsd.
* CMatrix.cc (ComplexMatrix::lssolve): Likewise, for zgelsd.
* fCMatrix.cc (FloatComplexMatrix::lssolve): Likewise, for cgelsd.
* test/bug-46330.tst: Don't avoid test on Windows systems.
| author | Avinoam Kalma <a.kalma@gmail.com> |
|---|---|
| date | Thu, 06 Oct 2016 11:52:12 -0400 |
| parents | 8d3a2d1af389 |
| children | 51b395d24782 |
| files | liboctave/array/CMatrix.cc liboctave/array/dMatrix.cc liboctave/array/fCMatrix.cc liboctave/array/fMatrix.cc test/bug-46330.tst |
| diffstat | 5 files changed, 91 insertions(+), 48 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/array/CMatrix.cc Wed Oct 05 17:35:24 2016 +0200 +++ b/liboctave/array/CMatrix.cc Thu Oct 06 11:52:12 2016 -0400 @@ -2415,6 +2415,7 @@ double dminmn = static_cast<double> (minmn); double dsmlsizp1 = static_cast<double> (smlsiz+1); double tmp = octave::math::log2 (dminmn / dsmlsizp1); + double anorm = 0.0; octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) @@ -2473,18 +2474,29 @@ lwork = static_cast<octave_idx_type> (octave::math::real (work(0))); work.resize (dim_vector (lwork, 1)); - F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, F77_DBLE_CMPLX_ARG (tmp_data), m, - F77_DBLE_CMPLX_ARG (pretval), - maxmn, ps, rcon, rank, - F77_DBLE_CMPLX_ARG (work.fortran_vec ()), lwork, - prwork, piwork, info)); - - if (s.elem (0) == 0.0) - rcon = 0.0; + anorm = xnorm (*this, 1); + + if (octave::math::isinf (anorm) || octave::math::isnan (anorm)) + { + rcon = 0.0; + octave::warn_singular_matrix (); + retval = Matrix (n, m, 0.0); + } else - rcon = s.elem (minmn - 1) / s.elem (0); - - retval.resize (n, nrhs); + { + F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, F77_DBLE_CMPLX_ARG (tmp_data), + m, F77_DBLE_CMPLX_ARG (pretval), + maxmn, ps, rcon, rank, + F77_DBLE_CMPLX_ARG (work.fortran_vec ()), + lwork, prwork, piwork, info)); + + if (s.elem (0) == 0.0) + rcon = 0.0; + else + rcon = s.elem (minmn - 1) / s.elem (0); + + retval.resize (n, nrhs); + } } return retval;
--- a/liboctave/array/dMatrix.cc Wed Oct 05 17:35:24 2016 +0200 +++ b/liboctave/array/dMatrix.cc Thu Oct 06 11:52:12 2016 -0400 @@ -2076,6 +2076,7 @@ double dminmn = static_cast<double> (minmn); double dsmlsizp1 = static_cast<double> (smlsiz+1); double tmp = octave::math::log2 (dminmn / dsmlsizp1); + double anorm = 0.0; octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) @@ -2131,17 +2132,28 @@ lwork = static_cast<octave_idx_type> (work(0)); work.resize (dim_vector (lwork, 1)); - F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcon, rank, - work.fortran_vec (), lwork, - piwork, info)); - - if (s.elem (0) == 0.0) - rcon = 0.0; + anorm = xnorm (*this, 1); + + if (octave::math::isinf (anorm) || octave::math::isnan (anorm)) + { + rcon = 0.0; + octave::warn_singular_matrix (); + retval = Matrix (n, m, 0.0); + } else - rcon = s.elem (minmn - 1) / s.elem (0); - - retval.resize (n, nrhs); + { + F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, + maxmn, ps, rcon, rank, + work.fortran_vec (), lwork, + piwork, info)); + + if (s.elem (0) == 0.0) + rcon = 0.0; + else + rcon = s.elem (minmn - 1) / s.elem (0); + + retval.resize (n, nrhs); + } } return retval;
--- a/liboctave/array/fCMatrix.cc Wed Oct 05 17:35:24 2016 +0200 +++ b/liboctave/array/fCMatrix.cc Thu Oct 06 11:52:12 2016 -0400 @@ -2432,6 +2432,7 @@ float dminmn = static_cast<float> (minmn); float dsmlsizp1 = static_cast<float> (smlsiz+1); float tmp = octave::math::log2 (dminmn / dsmlsizp1); + float anorm = 0.0; octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) @@ -2490,18 +2491,29 @@ lwork = static_cast<octave_idx_type> (octave::math::real (work(0))); work.resize (dim_vector (lwork, 1)); - F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, F77_CMPLX_ARG (tmp_data), m, - F77_CMPLX_ARG (pretval), - maxmn, ps, rcon, rank, - F77_CMPLX_ARG (work.fortran_vec ()), lwork, - prwork, piwork, info)); - - if (s.elem (0) == 0.0) - rcon = 0.0; + anorm = xnorm (*this, 1); + + if (octave::math::isinf (anorm) || octave::math::isnan (anorm)) + { + rcon = 0.0; + octave::warn_singular_matrix (); + retval = Matrix (n, m, 0.0); + } else - rcon = s.elem (minmn - 1) / s.elem (0); - - retval.resize (n, nrhs); + { + F77_XFCN (cgelsd, CGELSD, (m, n, nrhs, F77_CMPLX_ARG (tmp_data), + m, F77_CMPLX_ARG (pretval), + maxmn, ps, rcon, rank, + F77_CMPLX_ARG (work.fortran_vec ()), + lwork, prwork, piwork, info)); + + if (s.elem (0) == 0.0) + rcon = 0.0; + else + rcon = s.elem (minmn - 1) / s.elem (0); + + retval.resize (n, nrhs); + } } return retval;
--- a/liboctave/array/fMatrix.cc Wed Oct 05 17:35:24 2016 +0200 +++ b/liboctave/array/fMatrix.cc Thu Oct 06 11:52:12 2016 -0400 @@ -2103,6 +2103,7 @@ float dminmn = static_cast<float> (minmn); float dsmlsizp1 = static_cast<float> (smlsiz+1); float tmp = octave::math::log2 (dminmn / dsmlsizp1); + float anorm = 0.0; octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) @@ -2158,17 +2159,28 @@ lwork = static_cast<octave_idx_type> (work(0)); work.resize (dim_vector (lwork, 1)); - F77_XFCN (sgelsd, SGELSD, (m, n, nrhs, tmp_data, m, pretval, - maxmn, ps, rcon, rank, - work.fortran_vec (), lwork, - piwork, info)); - - if (s.elem (0) == 0.0) - rcon = 0.0; + anorm = xnorm (*this, 1); + + if (octave::math::isinf (anorm) || octave::math::isnan (anorm)) + { + rcon = 0.0; + octave::warn_singular_matrix (); + retval = Matrix (n, m, 0.0); + } else - rcon = s.elem (minmn - 1) / s.elem (0); - - retval.resize (n, nrhs); + { + F77_XFCN (sgelsd, SGELSD, (m, n, nrhs, tmp_data, m, pretval, + maxmn, ps, rcon, rank, + work.fortran_vec (), lwork, + piwork, info)); + + if (s.elem (0) == 0.0) + rcon = 0.0; + else + rcon = s.elem (minmn - 1) / s.elem (0); + + retval.resize (n, nrhs); + } } return retval;
--- a/test/bug-46330.tst Wed Oct 05 17:35:24 2016 +0200 +++ b/test/bug-46330.tst Thu Oct 06 11:52:12 2016 -0400 @@ -17,9 +17,4 @@ %! ## 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 +%! assert (c / (i * Inf * eye (4) - a) * b, zeros (2, 2))