Mercurial > hg > octave-nkf

changeset 15887:8ced82e96b48 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix segfaults with gesdd driver for svd (bug #37998). * liboctave/CmplxSVD.cc(init): Correctly size rwork array for gesdd driver. * liboctave/fCmplxSVD.cc(init): Correctly size rwork array for gesdd driver. * liboctave/dbleSVD.cc(init): Tweak coding style to match CmplxSVD.cc. * liboctave/floatSVD.cc(init): Tweak coding style to match fCmplxSVD.cc. * src/DLD-FUNCTIONS/svd.cc: Add %!test for gesdd driver and complex matrices.
author Rik <rik@octave.org>
date Thu, 03 Jan 2013 10:05:03 -0800
parents 065bc7944335
children fde9297ae074 7f4082b00a99
files liboctave/CmplxSVD.cc liboctave/dbleSVD.cc liboctave/fCmplxSVD.cc liboctave/floatSVD.cc src/DLD-FUNCTIONS/svd.cc
diffstat 5 files changed, 51 insertions(+), 24 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CmplxSVD.cc
+++ b/liboctave/CmplxSVD.cc
@@ -118,8 +118,7 @@
       //
       // For Lapack 3.0, this problem seems to be fixed.
 
-      jobu = 'N';
-      jobv = 'N';
+      jobu = jobv = 'N';
       ncol_u = nrow_vt = 1;
       break;
 
@@ -142,21 +141,21 @@
 
   Complex *vt = right_sm.fortran_vec ();
 
-  octave_idx_type lrwork = 5*max_mn;
-
-  Array<double> rwork (dim_vector (lrwork, 1));
-
-  // Ask ZGESVD what the dimension of WORK should be.
+  // Query ZGESVD for the correct dimension of WORK.
 
   octave_idx_type lwork = -1;
 
   Array<Complex> work (dim_vector (1, 1));
 
   octave_idx_type one = 1;
-  octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one);
+  octave_idx_type m1 = std::max (m, one);
+  octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
 
   if (svd_driver == SVD::GESVD)
     {
+      octave_idx_type lrwork = 5*max_mn;
+      Array<double> rwork (dim_vector (lrwork, 1));
+
       F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
                                  F77_CONST_CHAR_ARG2 (&jobv, 1),
                                  m, n, tmp_data, m1, s_vec, u, m1, vt,
@@ -180,6 +179,14 @@
     {
       assert (jobu == jobv);
       char jobz = jobu;
+
+      octave_idx_type lrwork;
+      if (jobz == 'N')
+        lrwork = 7*min_mn;
+      else
+        lrwork = 5*min_mn*min_mn + 5*min_mn;
+      Array<double> rwork (dim_vector (lrwork, 1));
+
       OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
 
       F77_XFCN (zgesdd, ZGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
--- a/liboctave/dbleSVD.cc
+++ b/liboctave/dbleSVD.cc
@@ -118,8 +118,7 @@
       //
       // For Lapack 3.0, this problem seems to be fixed.
 
-      jobu = 'N';
-      jobv = 'N';
+      jobu = jobv = 'N';
       ncol_u = nrow_vt = 1;
       break;
 
@@ -142,14 +141,15 @@
 
   double *vt = right_sm.fortran_vec ();
 
-  // Ask DGESVD what the dimension of WORK should be.
+  // Query DGESVD for the correct dimension of WORK.
 
   octave_idx_type lwork = -1;
 
   Array<double> work (dim_vector (1, 1));
 
   octave_idx_type one = 1;
-  octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one);
+  octave_idx_type m1 = std::max (m, one);
+  octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
 
   if (svd_driver == SVD::GESVD)
     {
--- a/liboctave/fCmplxSVD.cc
+++ b/liboctave/fCmplxSVD.cc
@@ -120,8 +120,7 @@
       //
       // For Lapack 3.0, this problem seems to be fixed.
 
-      jobu = 'N';
-      jobv = 'N';
+      jobu = jobv = 'N';
       ncol_u = nrow_vt = 1;
       break;
 
@@ -144,21 +143,21 @@
 
   FloatComplex *vt = right_sm.fortran_vec ();
 
-  octave_idx_type lrwork = 5*max_mn;
-
-  Array<float> rwork (dim_vector (lrwork, 1));
-
-  // Ask ZGESVD what the dimension of WORK should be.
+  // Query CGESVD for the correct dimension of WORK.
 
   octave_idx_type lwork = -1;
 
   Array<FloatComplex> work (dim_vector (1, 1));
 
   octave_idx_type one = 1;
-  octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one);
+  octave_idx_type m1 = std::max (m, one);
+  octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
 
   if (svd_driver == SVD::GESVD)
     {
+      octave_idx_type lrwork = 5*max_mn;
+      Array<float> rwork (dim_vector (lrwork, 1));
+
       F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
                                  F77_CONST_CHAR_ARG2 (&jobv, 1),
                                  m, n, tmp_data, m1, s_vec, u, m1, vt,
@@ -182,6 +181,14 @@
     {
       assert (jobu == jobv);
       char jobz = jobu;
+
+      octave_idx_type lrwork;
+      if (jobz == 'N')
+        lrwork = 5*min_mn;
+      else
+        lrwork = min_mn * std::max (5*min_mn+7, 2*max_mn+2*min_mn+1);
+      Array<float> rwork (dim_vector (lrwork, 1));
+
       OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
 
       F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
--- a/liboctave/floatSVD.cc
+++ b/liboctave/floatSVD.cc
@@ -118,8 +118,7 @@
       //
       // For Lapack 3.0, this problem seems to be fixed.
 
-      jobu = 'N';
-      jobv = 'N';
+      jobu = jobv = 'N';
       ncol_u = nrow_vt = 1;
       break;
 
@@ -142,14 +141,15 @@
 
   float *vt = right_sm.fortran_vec ();
 
-  // Ask DGESVD what the dimension of WORK should be.
+  // Query SGESVD for the correct dimension of WORK.
 
   octave_idx_type lwork = -1;
 
   Array<float> work (dim_vector (1, 1));
 
   octave_idx_type one = 1;
-  octave_idx_type m1 = std::max (m, one), nrow_vt1 = std::max (nrow_vt, one);
+  octave_idx_type m1 = std::max (m, one);
+  octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
 
   if (svd_driver == SVD::GESVD)
     {
--- a/src/DLD-FUNCTIONS/svd.cc
+++ b/src/DLD-FUNCTIONS/svd.cc
@@ -423,3 +423,16 @@
 
   return SET_INTERNAL_VARIABLE_CHOICES (svd_driver, driver_names);
 }
+
+/*
+%!test
+%! A = [1+1i, 1-1i, 0; 0, 2, 0; 1i, 1i, 1+2i];
+%! old_driver = svd_driver ("gesvd");
+%! [U1, S1, V1] = svd (A);
+%! svd_driver ("gesdd");
+%! [U2, S2, V2] = svd (A);
+%! assert (U1, U2, 5*eps);
+%! assert (S1, S2, 5*eps);
+%! assert (V1, V2, 5*eps);
+%! svd_driver (old_driver);
+*/