# HG changeset patch # User Rik # Date 1357236303 28800 # Node ID 8ced82e96b48d9d98dcce95295c8d901feac4090 # Parent 065bc79443350b77397102c9febb5dab85a61c1a 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. diff -r 065bc7944335 -r 8ced82e96b48 liboctave/CmplxSVD.cc --- a/liboctave/CmplxSVD.cc Thu Jan 03 11:16:40 2013 -0500 +++ b/liboctave/CmplxSVD.cc Thu Jan 03 10:05:03 2013 -0800 @@ -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 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 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 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 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), diff -r 065bc7944335 -r 8ced82e96b48 liboctave/dbleSVD.cc --- a/liboctave/dbleSVD.cc Thu Jan 03 11:16:40 2013 -0500 +++ b/liboctave/dbleSVD.cc Thu Jan 03 10:05:03 2013 -0800 @@ -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 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) { diff -r 065bc7944335 -r 8ced82e96b48 liboctave/fCmplxSVD.cc --- a/liboctave/fCmplxSVD.cc Thu Jan 03 11:16:40 2013 -0500 +++ b/liboctave/fCmplxSVD.cc Thu Jan 03 10:05:03 2013 -0800 @@ -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 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 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 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 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), diff -r 065bc7944335 -r 8ced82e96b48 liboctave/floatSVD.cc --- a/liboctave/floatSVD.cc Thu Jan 03 11:16:40 2013 -0500 +++ b/liboctave/floatSVD.cc Thu Jan 03 10:05:03 2013 -0800 @@ -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 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) { diff -r 065bc7944335 -r 8ced82e96b48 src/DLD-FUNCTIONS/svd.cc --- a/src/DLD-FUNCTIONS/svd.cc Thu Jan 03 11:16:40 2013 -0500 +++ b/src/DLD-FUNCTIONS/svd.cc Thu Jan 03 10:05:03 2013 -0800 @@ -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); +*/