Mercurial > octave
changeset 27093:6e18f0ce268c
Inverse of a sparse/diagonal singular matrix should be a sparse/diagonal matrix of Infs.
* liboctave/array/CSparse.cc(inverse): if the matrix is null, return an error.
If it is singular, return a sparse matrix of Infs.
* liboctave/array/dSparse.cc(inverse): if the matrix is null, return an error.
If it is singular, return a sparse matrix of Infs.
* liboctave/array/CDiagMatrix.cc(inverse): if the matrix is null, return
an error. If it is singular, return a diagonal matrix of Infs.
* liboctave/array/dDiagMatrix.cc(inverse): if the matrix is null, return
an error. If it is singular, return a diagonal matrix of Infs.
* libinterp/corefcn/inv.cc: set rcond to zero for singular diagonal matrices.
Add tests for the above cases.
* test/diag-perm.tst: modify a test for the inversion of a singular diagonal
matrix.
| author | marco.caliari@univr.it |
|---|---|
| date | Tue, 14 May 2019 17:00:52 +0200 |
| parents | d1fe5bd5e005 |
| children | f16471efcdf4 |
| files | libinterp/corefcn/inv.cc liboctave/array/CDiagMatrix.cc liboctave/array/CSparse.cc liboctave/array/dDiagMatrix.cc liboctave/array/dSparse.cc test/diag-perm.tst |
| diffstat | 6 files changed, 124 insertions(+), 17 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/inv.cc Sat May 18 13:18:00 2019 +0200 +++ b/libinterp/corefcn/inv.cc Tue May 14 17:00:52 2019 +0200 @@ -83,13 +83,17 @@ if (isfloat) { result = arg.float_complex_diag_matrix_value ().inverse (info); - if (nargout > 1) + if (info == -1) + frcond = 0.0f; + else if (nargout > 1) frcond = arg.float_complex_diag_matrix_value ().rcond (); } else { result = arg.complex_diag_matrix_value ().inverse (info); - if (nargout > 1) + if (info == -1) + rcond = 0.0; + else if (nargout > 1) rcond = arg.complex_diag_matrix_value ().rcond (); } } @@ -98,13 +102,17 @@ if (isfloat) { result = arg.float_diag_matrix_value ().inverse (info); - if (nargout > 1) + if (info == -1) + frcond = 0.0f; + else if (nargout > 1) frcond = arg.float_diag_matrix_value ().rcond (); } else { result = arg.diag_matrix_value ().inverse (info); - if (nargout > 1) + if (info == -1) + rcond = 0.0; + else if (nargout > 1) rcond = arg.diag_matrix_value ().rcond (); } } @@ -189,7 +197,7 @@ if (isfloat) { volatile float xrcond = frcond; - rcond_plus_one_eq_one = xrcond + 1.0F == 1.0F; + rcond_plus_one_eq_one = xrcond + 1.0f == 1.0f; } else { @@ -225,10 +233,29 @@ %! assert (isa (xinv, "double")); %! assert (isa (rcond, "double")); +%!test +%! fail ("A = inv (sparse ([1, 2;0 ,0]))", "warning", "matrix singular"); +%! assert (A, sparse ([Inf, Inf; 0, 0])); + +%!test +%! 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]))) + +%!test +%! 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])); + %!error inv () %!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)) */ DEFALIAS (inverse, inv);
--- a/liboctave/array/CDiagMatrix.cc Sat May 18 13:18:00 2019 +0200 +++ b/liboctave/array/CDiagMatrix.cc Tue May 14 17:00:52 2019 +0200 @@ -316,15 +316,35 @@ ComplexDiagMatrix retval (r, c); info = 0; - for (octave_idx_type i = 0; i < length (); i++) + octave_idx_type len = r; // alias for readability + octave_idx_type z_count = 0; // zeros + octave_idx_type nz_count = 0; // non-zeros + for (octave_idx_type i = 0; i < len; i++) { - if (elem (i, i) == 0.0) + if (xelem (i, i) == 0.0) { - info = -1; - return *this; + z_count++; + if (nz_count > 0) + break; } else - retval.elem (i, i) = 1.0 / elem (i, i); + { + nz_count++; + if (z_count > 0) + break; + retval.elem (i, i) = 1.0 / xelem (i, i); + } + } + if (nz_count == 0) + { + (*current_liboctave_error_handler) + ("inverse of the null matrix not defined"); + } + else if (z_count > 0) + { + info = -1; + element_type *data = retval.fortran_vec (); + std::fill (data, data + len, octave::numeric_limits<double>::Inf ()); } return retval;
--- a/liboctave/array/CSparse.cc Sat May 18 13:18:00 2019 +0200 +++ b/liboctave/array/CSparse.cc Tue May 14 17:00:52 2019 +0200 @@ -985,6 +985,12 @@ SparseComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info, double& rcond, bool, bool calc_cond) const { + if (nnz () == 0) + { + (*current_liboctave_error_handler) + ("inverse of the null matrix not defined"); + } + int typ = mattype.type (false); SparseComplexMatrix ret; @@ -1035,6 +1041,18 @@ Qinit, Matrix (), false, false); rcond = fact.rcond (); + if (rcond == 0.0) + { + // Return all Inf matrix with sparsity pattern of input. + octave_idx_type nz = nnz (); + ret = SparseComplexMatrix (rows (), cols (), nz); + std::fill (ret.xdata (), ret.xdata () + nz, + octave::numeric_limits<double>::Inf ()); + std::copy_n (ridx (), nz, ret.xridx ()); + std::copy_n (cidx (), cols () + 1, ret.xcidx ()); + + return ret; + } double rcond2; SparseComplexMatrix InvL = fact.L ().transpose (). tinverse (tmp_typ, info, rcond2,
--- a/liboctave/array/dDiagMatrix.cc Sat May 18 13:18:00 2019 +0200 +++ b/liboctave/array/dDiagMatrix.cc Tue May 14 17:00:52 2019 +0200 @@ -233,19 +233,41 @@ { octave_idx_type r = rows (); octave_idx_type c = cols (); - octave_idx_type len = length (); if (r != c) (*current_liboctave_error_handler) ("inverse requires square matrix"); DiagMatrix retval (r, c); info = 0; + octave_idx_type len = r; // alias for readability + octave_idx_type z_count = 0; // zeros + octave_idx_type nz_count = 0; // non-zeros for (octave_idx_type i = 0; i < len; i++) { - if (elem (i, i) == 0.0) - retval.elem (i, i) = octave::numeric_limits<double>::Inf (); + if (xelem (i, i) == 0.0) + { + z_count++; + if (nz_count > 0) + break; + } else - retval.elem (i, i) = 1.0 / elem (i, i); + { + nz_count++; + if (z_count > 0) + break; + retval.elem (i, i) = 1.0 / xelem (i, i); + } + } + if (nz_count == 0) + { + (*current_liboctave_error_handler) + ("inverse of the null matrix not defined"); + } + else if (z_count > 0) + { + info = -1; + element_type *data = retval.fortran_vec (); + std::fill (data, data + len, octave::numeric_limits<double>::Inf ()); } return retval;
--- a/liboctave/array/dSparse.cc Sat May 18 13:18:00 2019 +0200 +++ b/liboctave/array/dSparse.cc Tue May 14 17:00:52 2019 +0200 @@ -928,6 +928,12 @@ SparseMatrix::inverse (MatrixType& mattype, octave_idx_type& info, double& rcond, bool, bool calc_cond) const { + if (nnz () == 0) + { + (*current_liboctave_error_handler) + ("inverse of the null matrix not defined"); + } + int typ = mattype.type (false); SparseMatrix ret; @@ -977,6 +983,19 @@ Qinit, Matrix (), false, false); rcond = fact.rcond (); + if (rcond == 0.0) + { + // Return all Inf matrix with sparsity pattern of input. + octave_idx_type nz = nnz (); + ret = SparseMatrix (rows (), cols (), nz); + std::fill (ret.xdata (), ret.xdata () + nz, + octave::numeric_limits<double>::Inf ()); + std::copy_n (ridx (), nz, ret.xridx ()); + std::copy_n (cidx (), cols () + 1, ret.xcidx ()); + + return ret; + } + double rcond2; SparseMatrix InvL = fact.L ().transpose ().tinverse (tmp_typ, info, rcond2, true, false);
--- a/test/diag-perm.tst Sat May 18 13:18:00 2019 +0200 +++ b/test/diag-perm.tst Tue May 14 17:00:52 2019 +0200 @@ -265,11 +265,12 @@ %! assert (full (D - A), D - full (A)); ## inverse preserves diagonal structure even for singular matrices (bug #46103) +## but set all the diagonal elements to Inf (bug #56232) %!test %! x = diag (1:3); %! assert (inv (x), diag ([1 1/2 1/3])); -%! x = diag (0:2); -%! assert (inv (x), diag ([Inf 1 1/2])); +%!warning <matrix singular> A = inv (diag (0:2)); +%! assert (A, diag ([Inf Inf Inf])); ## assignment to diagonal elements preserves diagonal structure (bug #36932) %!test