Mercurial > hg > octave-nkf
changeset 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 | 9c782a34b545 |
files | liboctave/ChangeLog liboctave/CmplxSVD.cc liboctave/dbleSVD.cc |
diffstat | 3 files changed, 50 insertions(+), 21 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,5 +1,8 @@ 1999-11-03 John W. Eaton <jwe@bevo.che.wisc.edu> + * dbleSVD.cc (SVD::init): Let DGESVD determine work space requirement. + * CmplxSVD.cc (ComplexSVD::init): Likewise, for complex version. + * dbleSCHUR.cc (SCHUR::init): IWORK is always referenced by dgeesx. Don't forget to pass length of third string argument to dgeesx.
--- a/liboctave/CmplxSVD.cc +++ b/liboctave/CmplxSVD.cc @@ -131,26 +131,41 @@ Complex *vt = right_sm.fortran_vec (); - int lwork = 2*min_mn + max_mn; - - Array<Complex> work (lwork); - Complex *pwork = work.fortran_vec (); - int lrwork = 5*max_mn; Array<double> rwork (lrwork); - double *prwork = rwork.fortran_vec (); + + // Ask ZGESVD what the dimension of WORK should be. + + int lwork = -1; - F77_XFCN (zgesvd, ZGESVD, (&jobu, &jobv, m, n, tmp_data, m, s_vec, u, - m, vt, nrow_vt, pwork, lwork, prwork, info, - 1L, 1L)); + Array<Complex> work (1); + + F77_XFCN (zgesvd, ZGESVD, (&jobu, &jobv, m, n, tmp_data, m, s_vec, + u, m, vt, nrow_vt, work.fortran_vec (), + lwork, rwork.fortran_vec (), info, 1L, + 1L)); if (f77_exception_encountered) (*current_liboctave_error_handler) ("unrecoverable error in zgesvd"); else { - if (! (jobv == 'N' || jobv == 'O')) - right_sm = right_sm.hermitian (); + lwork = static_cast<int> (work(0).real ()); + work.resize (lwork); + + F77_XFCN (zgesvd, ZGESVD, (&jobu, &jobv, m, n, tmp_data, m, + s_vec, u, m, vt, nrow_vt, + work.fortran_vec (), lwork, + rwork.fortran_vec (), + info, 1L, 1L)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) ("unrecoverable error in zgesvd"); + else + { + if (! (jobv == 'N' || jobv == 'O')) + right_sm = right_sm.hermitian (); + } } return info;
--- 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;