Mercurial > hg > octave-nkf

changeset 10185:455759a5fcbe

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 ())