From owner-bug-octave at bevo dot che dot wisc dot edu Wed Jan 22 16:33:49 1997 Subject: ovtave-2.0.1: trouble with svd From: "John W. Eaton" To: pfost at aghl dot ruhr-uni-bochum dot de (Martin Pfost) Cc: bug-octave at bevo dot che dot wisc dot edu Date: Wed, 22 Jan 1997 16:32:42 -0600 On 14-Jan-1997, Martin Pfost wrote: : it seems to me that the singular value decomposition seems is buggy on : both Linux and HP-UX. As far as I can tell, the following result : should be non-negative: : : octave:1> svd([eye(3) eye(3)]) : ans = : : -1.4142 : -1.4142 : -1.4142 : : This also messes up the rank()-function: : : octave:3> rank([eye(3) eye(3)]) : ans = 0 : : However, if called with three return values, I get the correct results: This seems to be a bug in Lapack. It appers in both dgesvd() and zgesvd(). I suspect it also appears in sgesvd() and cgesvd(), but I haven't actually tested them. Here is a simple test program that demonstrates the problem. program main integer info double precision a(3,6), s(3), u(3,3), vt(6,6), work(1000) data a /18*0.0/ a(1,1) = 1.0 a(2,2) = 1.0 a(3,3) = 1.0 a(1,4) = 1.0 a(2,5) = 1.0 a(3,6) = 1.0 * This call should return s = [1.41421356 1.41421356 1.41421356] * but with Lapack 2.0, it returns [-1.41421356 -1.41421356 -1.41421356]. call dgesvd ('N', 'N', 3, 6, a, 3, s, u, 3, vt, 6, $ work, 1000, info) print *, 's' print *, s * This one does return s = [1.41421356 1.41421356 1.41421356]. call dgesvd ('O', 'N', 3, 6, a, 3, s, u, 3, vt, 6, $ work, 1000, info) print *, 's' print *, s end Unless someone can show me that there is a bug in the way I am calling dgesvd for the case of jobu == jobvt == 'N', then I have to conclude that the bug is in Lapack. In any case, the following patch for Octave appears to work around the problem. Thanks jwe Wed Jan 22 16:18:53 1997 John W. Eaton * dbleSVD.cc (SVD::init): Work around apparent dgesvd() bug. * CmplxSVD.cc (ComplexSVD::init): Work around apparent zgesvd() bug. Index: CmplxSVD.cc =================================================================== RCS file: /home/jwe/src/master/octave/liboctave/CmplxSVD.cc,v retrieving revision 1.18 diff -c -r1.18 CmplxSVD.cc *** CmplxSVD.cc 1996/03/03 01:16:15 1.18 --- CmplxSVD.cc 1997/01/22 22:11:09 *************** *** 99,105 **** break; case SVD::sigma_only: ! jobu = jobv = 'N'; ncol_u = nrow_vt = 1; break; --- 99,113 ---- break; case SVD::sigma_only: ! ! // Note: for this case, both jobu and jobv should be 'N', but ! // there seems to be a bug in dgesvd from Lapack V2.0. To ! // demonstrate the bug, set both jobu and jobv to 'N' and find ! // the singular values of [eye(3), eye(3)]. The result is ! // [-sqrt(2), -sqrt(2), -sqrt(2)]. ! ! jobu = 'O'; ! jobv = 'N'; ncol_u = nrow_vt = 1; break; *************** *** 109,115 **** type_computed = svd_type; ! if (jobu != 'N') left_sm.resize (m, ncol_u); Complex *u = left_sm.fortran_vec (); --- 117,123 ---- type_computed = svd_type; ! if (! (jobu == 'N' || jobu == 'O')) left_sm.resize (m, ncol_u); Complex *u = left_sm.fortran_vec (); *************** *** 117,123 **** sigma.resize (nrow_s, ncol_s); double *s_vec = sigma.fortran_vec (); ! if (jobv != 'N') right_sm.resize (nrow_vt, n); Complex *vt = right_sm.fortran_vec (); --- 125,131 ---- sigma.resize (nrow_s, ncol_s); double *s_vec = sigma.fortran_vec (); ! if (! (jobv == 'N' || jobv == 'O')) right_sm.resize (nrow_vt, n); Complex *vt = right_sm.fortran_vec (); *************** *** 140,146 **** (*current_liboctave_error_handler) ("unrecoverable error in zgesvd"); else { ! if (jobv != 'N') right_sm = right_sm.hermitian (); } --- 148,154 ---- (*current_liboctave_error_handler) ("unrecoverable error in zgesvd"); else { ! if (! (jobv == 'N' || jobv == 'O')) right_sm = right_sm.hermitian (); } Index: dbleSVD.cc =================================================================== RCS file: /home/jwe/src/master/octave/liboctave/dbleSVD.cc,v retrieving revision 1.19 diff -c -r1.19 dbleSVD.cc *** dbleSVD.cc 1996/03/03 01:16:15 1.19 --- dbleSVD.cc 1997/01/22 22:11:14 *************** *** 99,105 **** break; case SVD::sigma_only: ! jobu = jobv = 'N'; ncol_u = nrow_vt = 1; break; --- 99,113 ---- break; case SVD::sigma_only: ! ! // Note: for this case, both jobu and jobv should be 'N', but ! // there seems to be a bug in dgesvd from Lapack V2.0. To ! // demonstrate the bug, set both jobu and jobv to 'N' and find ! // the singular values of [eye(3), eye(3)]. The result is ! // [-sqrt(2), -sqrt(2), -sqrt(2)]. ! ! jobu = 'O'; ! jobv = 'N'; ncol_u = nrow_vt = 1; break; *************** *** 109,115 **** type_computed = svd_type; ! if (jobu != 'N') left_sm.resize (m, ncol_u); double *u = left_sm.fortran_vec (); --- 117,123 ---- type_computed = svd_type; ! if (! (jobu == 'N' || jobu == 'O')) left_sm.resize (m, ncol_u); double *u = left_sm.fortran_vec (); *************** *** 117,123 **** sigma.resize (nrow_s, ncol_s); double *s_vec = sigma.fortran_vec (); ! if (jobv != 'N') right_sm.resize (nrow_vt, n); double *vt = right_sm.fortran_vec (); --- 125,131 ---- sigma.resize (nrow_s, ncol_s); double *s_vec = sigma.fortran_vec (); ! if (! (jobv == 'N' || jobv == 'O')) right_sm.resize (nrow_vt, n); double *vt = right_sm.fortran_vec (); *************** *** 137,143 **** (*current_liboctave_error_handler) ("unrecoverable error in dgesvd"); else { ! if (jobv != 'N') right_sm = right_sm.transpose (); } --- 145,151 ---- (*current_liboctave_error_handler) ("unrecoverable error in dgesvd"); else { ! if (! (jobv == 'N' || jobv == 'O')) right_sm = right_sm.transpose (); }