Mercurial > hg > octave-nkf
changeset 11668:8ac2994e4596 release-3-0-x
more xGELSD workspace fixes
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Mon, 03 Mar 2008 02:22:03 -0500 |
parents | bc4b8d973a3a |
children | a6e08ecb4050 |
files | liboctave/CMatrix.cc liboctave/ChangeLog liboctave/dMatrix.cc |
diffstat | 3 files changed, 31 insertions(+), 16 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -2520,11 +2520,11 @@ double dminmn = static_cast<double> (minmn); double dsmlsizp1 = static_cast<double> (smlsiz+1); #if defined (HAVE_LOG2) - double tmp = log2 (dminmn) / dsmlsizp1 + 1; + double tmp = log2 (dminmn / dsmlsizp1); #else - double tmp = log (dminmn) / dsmlsizp1 / log (2.0) + 1; + double tmp = log (dminmn / dsmlsizp1) / log (2.0); #endif - octave_idx_type nlvl = static_cast<int> (tmp); + octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) nlvl = 0; @@ -2713,8 +2713,12 @@ Array<Complex> work (1); - // FIXME: Can SMLSIZ be other than 25? - octave_idx_type smlsiz = 25; + octave_idx_type smlsiz; + F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("ZGELSD", 6), + F77_CONST_CHAR_ARG2 (" ", 1), + 0, 0, 0, 0, smlsiz + F77_CHAR_ARG_LEN (6) + F77_CHAR_ARG_LEN (1)); // We compute the size of rwork and iwork because ZGELSD in // older versions of LAPACK does not return them on a query @@ -2722,11 +2726,11 @@ double dminmn = static_cast<double> (minmn); double dsmlsizp1 = static_cast<double> (smlsiz+1); #if defined (HAVE_LOG2) - double tmp = log2 (dminmn) / dsmlsizp1 + 1; + double tmp = log2 (dminmn / dsmlsizp1); #else - double tmp = log (dminmn) / dsmlsizp1 / log (2.0) + 1; + double tmp = log (dminmn / dsmlsizp1) / log (2.0); #endif - octave_idx_type nlvl = static_cast<int> (tmp); + octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) nlvl = 0;
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,10 @@ +2008-03-03 Jaroslav Hajek <highegg@gmail.com> + + * dMatrix.cc (Matrix::lssolve): Also avoid dgelsd lwork query bug + in lssolve method that accepts column vector argument. Correct + calculation of nlvl. + * CMatrix.cc (ComplexMatrix::lssolve): Likewise, for zgelsd. + 2008-02-26 John W. Eaton <jwe@octave.org> * oct-rand.cc (get_dist_id): Fix typo.
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -2130,11 +2130,11 @@ double dminmn = static_cast<double> (minmn); double dsmlsizp1 = static_cast<double> (smlsiz+1); #if defined (HAVE_LOG2) - double tmp = log2 (dminmn) / dsmlsizp1 + 1; + double tmp = log2 (dminmn) / dsmlsizp1; #else - double tmp = log (dminmn) / dsmlsizp1 / log (2.0) + 1; + double tmp = log (dminmn) / dsmlsizp1 / log (2.0); #endif - octave_idx_type nlvl = static_cast<int> (tmp); + octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) nlvl = 0; @@ -2324,19 +2324,23 @@ Array<double> work (1); - // FIXME: Can SMLSIZ be other than 25? - octave_idx_type smlsiz = 25; + octave_idx_type smlsiz; + F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6), + F77_CONST_CHAR_ARG2 (" ", 1), + 0, 0, 0, 0, smlsiz + F77_CHAR_ARG_LEN (6) + F77_CHAR_ARG_LEN (1)); // We compute the size of iwork because DGELSD in older versions // of LAPACK does not return it on a query call. double dminmn = static_cast<double> (minmn); double dsmlsizp1 = static_cast<double> (smlsiz+1); #if defined (HAVE_LOG2) - double tmp = log2 (dminmn) / dsmlsizp1 + 1; + double tmp = log2 (dminmn / dsmlsizp1); #else - double tmp = log (dminmn) / dsmlsizp1 / log (2.0) + 1; + double tmp = log (dminmn / dsmlsizp1) / log (2.0); #endif - octave_idx_type nlvl = static_cast<int> (tmp); + octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1; if (nlvl < 0) nlvl = 0;