Mercurial > hg > octave-lyh
diff liboctave/dbleSVD.cc @ 3336:08ad797989f8
[project @ 1999-11-03 21:41:34 by jwe]
author | jwe |
---|---|
date | Wed, 03 Nov 1999 21:41:35 +0000 |
parents | f39b97e13cf2 |
children | d14c483b3c12 |
line wrap: on
line diff
--- a/liboctave/dbleSVD.cc +++ b/liboctave/dbleSVD.cc @@ -80,7 +80,6 @@ double *tmp_data = atmp.fortran_vec (); int min_mn = m < n ? m : n; - int max_mn = m > n ? m : n; char jobu = 'A'; char jobv = 'A'; @@ -131,23 +130,35 @@ double *vt = right_sm.fortran_vec (); - int tmp1 = 3*min_mn + max_mn; - int tmp2 = 5*min_mn - 4; - int lwork = tmp1 > tmp2 ? tmp1 : tmp2; + // Ask DGESVD what the dimension of WORK should be. - Array<double> work (lwork); - double *pwork = work.fortran_vec (); + int lwork = -1; + + Array<double> work (1); F77_XFCN (dgesvd, DGESVD, (&jobu, &jobv, m, n, tmp_data, m, s_vec, - u, m, vt, nrow_vt, pwork, lwork, info, - 1L, 1L)); + u, m, vt, nrow_vt, work.fortran_vec (), + lwork, info, 1L, 1L)); if (f77_exception_encountered) (*current_liboctave_error_handler) ("unrecoverable error in dgesvd"); else { - if (! (jobv == 'N' || jobv == 'O')) - right_sm = right_sm.transpose (); + lwork = static_cast<int> (work(0)); + work.resize (lwork); + + F77_XFCN (dgesvd, DGESVD, (&jobu, &jobv, m, n, tmp_data, m, + s_vec, u, m, vt, nrow_vt, + work.fortran_vec (), lwork, info, 1L, + 1L)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in dgesvd"); + else + { + if (! (jobv == 'N' || jobv == 'O')) + right_sm = right_sm.transpose (); + } } return info;