Mercurial > hg > octave-terminal
changeset 10185:455759a5fcbe
fix norm and svd on empty matrices
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Fri, 22 Jan 2010 12:37:03 +0100 |
parents | b39bd23019eb |
children | 095a1e670e68 |
files | liboctave/ChangeLog liboctave/CmplxSVD.cc liboctave/dbleSVD.cc liboctave/fCmplxSVD.cc liboctave/floatSVD.cc src/ChangeLog src/DLD-FUNCTIONS/svd.cc src/data.cc |
diffstat | 8 files changed, 89 insertions(+), 40 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,10 @@ +2010-01-22 Jaroslav Hajek <highegg@gmail.com> + + * dbleSVD.cc (SVD::init): Ensure LD* arguments >= 1. + * floatSVD.cc (FloatSVD::init): Likewise. + * CmplxSVD.cc (ComplexSVD::init): Likewise. + * fCmplxSVD.cc (FloatComplexSVD::init): Likewise. + 2010-01-21 John W. Eaton <jwe@octave.org> * CMatrix.cc, fCMatrix.cc, cmd-edit.cc, cmd-hist.cc,
--- a/liboctave/CmplxSVD.cc +++ b/liboctave/CmplxSVD.cc @@ -141,10 +141,12 @@ Array<Complex> work (1); + octave_idx_type m1 = std::max (m, 1), nrow_vt1 = std::max (nrow_vt, 1); + F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m, s_vec, u, m, vt, - nrow_vt, work.fortran_vec (), lwork, + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -154,8 +156,8 @@ F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m, s_vec, u, m, vt, - nrow_vt, work.fortran_vec (), lwork, + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)));
--- a/liboctave/dbleSVD.cc +++ b/liboctave/dbleSVD.cc @@ -137,10 +137,12 @@ Array<double> work (1); + octave_idx_type m1 = std::max (m, 1), nrow_vt1 = std::max (nrow_vt, 1); + F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m, s_vec, u, m, vt, - nrow_vt, work.fortran_vec (), lwork, info + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -149,8 +151,8 @@ F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m, s_vec, u, m, vt, - nrow_vt, work.fortran_vec (), lwork, info + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)));
--- a/liboctave/fCmplxSVD.cc +++ b/liboctave/fCmplxSVD.cc @@ -141,10 +141,12 @@ Array<FloatComplex> work (1); + octave_idx_type m1 = std::max (m, 1), nrow_vt1 = std::max (nrow_vt, 1); + F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m, s_vec, u, m, vt, - nrow_vt, work.fortran_vec (), lwork, + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -154,8 +156,8 @@ F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m, s_vec, u, m, vt, - nrow_vt, work.fortran_vec (), lwork, + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, rwork.fortran_vec (), info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)));
--- a/liboctave/floatSVD.cc +++ b/liboctave/floatSVD.cc @@ -137,10 +137,12 @@ Array<float> work (1); + octave_idx_type m1 = std::max (m, 1), nrow_vt1 = std::max (nrow_vt, 1); + F77_XFCN (sgesvd, SGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m, s_vec, u, m, vt, - nrow_vt, work.fortran_vec (), lwork, info + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1))); @@ -149,8 +151,8 @@ F77_XFCN (sgesvd, SGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1), F77_CONST_CHAR_ARG2 (&jobv, 1), - m, n, tmp_data, m, s_vec, u, m, vt, - nrow_vt, work.fortran_vec (), lwork, info + m, n, tmp_data, m1, s_vec, u, m1, vt, + nrow_vt1, work.fortran_vec (), lwork, info F77_CHAR_ARG_LEN (1) F77_CHAR_ARG_LEN (1)));
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,8 @@ +2010-01-22 Jaroslav Hajek <highegg@gmail.com> + + * data.cc (Fnorm): Don't special-case empty matrix. + * DLD-FUNCTIONS/svd.cc (Fsvd): Fix empty matrix case. Add tests. + 2010-01-21 Jaroslav Hajek <highegg@gmail.com> * ov-builtin.cc (octave_builtin::do_multi_index_op): Partially undo
--- a/src/DLD-FUNCTIONS/svd.cc +++ b/src/DLD-FUNCTIONS/svd.cc @@ -134,39 +134,61 @@ octave_idx_type nr = arg.rows (); octave_idx_type nc = arg.columns (); + if (arg.ndims () != 2) + { + error ("svd: only valid for matrices"); + return retval; + } + bool isfloat = arg.is_single_type (); + SVD::type type = ((nargout == 0 || nargout == 1) + ? SVD::sigma_only + : (nargin == 2) ? SVD::economy : SVD::std); + if (nr == 0 || nc == 0) { if (isfloat) { - if (nargout == 3) + switch (type) { - retval(3) = float_identity_matrix (nr, nr); - retval(2) = FloatMatrix (nr, nc); - retval(1) = float_identity_matrix (nc, nc); + case SVD::std: + retval(2) = FloatDiagMatrix (nc, nc, 1.0f); + retval(1) = FloatMatrix (nr, nc); + retval(0) = FloatDiagMatrix (nr, nr, 1.0f); + break; + case SVD::economy: + retval(2) = FloatDiagMatrix (0, nc, 1.0f); + retval(1) = FloatMatrix (0, 0); + retval(0) = FloatDiagMatrix (nr, 0, 1.0f); + break; + case SVD::sigma_only: default: + retval(0) = FloatMatrix (0, 1); + break; } - else - retval(0) = FloatMatrix (0, 1); } else { - if (nargout == 3) + switch (type) { - retval(3) = identity_matrix (nr, nr); - retval(2) = Matrix (nr, nc); - retval(1) = identity_matrix (nc, nc); + case SVD::std: + retval(2) = DiagMatrix (nc, nc, 1.0); + retval(1) = Matrix (nr, nc); + retval(0) = DiagMatrix (nr, nr, 1.0); + break; + case SVD::economy: + retval(2) = DiagMatrix (0, nc, 1.0); + retval(1) = Matrix (0, 0); + retval(0) = DiagMatrix (nr, 0, 1.0); + break; + case SVD::sigma_only: default: + retval(0) = Matrix (0, 1); + break; } - else - retval(0) = Matrix (0, 1); } } else { - SVD::type type = ((nargout == 0 || nargout == 1) - ? SVD::sigma_only - : (nargin == 2) ? SVD::economy : SVD::std); - if (isfloat) { if (arg.is_real_type ()) @@ -355,6 +377,20 @@ %! [u, s, v] = svd (a, 1); %! assert (u * s * v', a, sqrt (eps('single'))); +%!test +%! a = zeros (0, 5); +%! [u, s, v] = svd (a); +%! assert (size (u), [0, 0]); +%! assert (size (s), [0, 5]); +%! assert (size (v), [5, 5]); + +%!test +%! a = zeros (5, 0); +%! [u, s, v] = svd (a, 1); +%! assert (size (u), [5, 0]); +%! assert (size (s), [0, 0]); +%! assert (size (v), [0, 0]); + %!error <Invalid call to svd.*> svd (); %!error <Invalid call to svd.*> svd ([1, 2; 4, 5], 2, 3); %!error <Invalid call to svd.*> [u, v] = svd ([1, 2; 3, 4]);
--- a/src/data.cc +++ b/src/data.cc @@ -5061,14 +5061,7 @@ { octave_value x_arg = args(0); - if (x_arg.is_empty ()) - { - if (x_arg.is_single_type ()) - retval(0) = static_cast<float>(0.0); - else - retval(0) = 0.0; - } - else if (x_arg.ndims () == 2) + if (x_arg.ndims () == 2) { enum { sfmatrix, sfcols, sfrows, sffrob, sfinf } strflag = sfmatrix; if (nargin > 1 && args(nargin-1).is_string ())