Mercurial > hg > octave-lyh
changeset 3752:719a44ff67c9
[project @ 2000-12-13 19:02:42 by jwe]
| author | jwe |
|---|---|
| date | Wed, 13 Dec 2000 19:02:43 +0000 |
| parents | 1ae5be669422 |
| children | f751e43de300 |
| files | liboctave/CMatrix.cc liboctave/ChangeLog liboctave/dMatrix.cc |
| diffstat | 3 files changed, 107 insertions(+), 72 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -1531,34 +1531,44 @@ double rcond = -1.0; - int lwork; - if (m < n) - lwork = 2*m + (nrhs > n ? nrhs : n); - else - lwork = 2*n + (nrhs > m ? nrhs : m); - - lwork *= 16; - - Array<Complex> work (lwork); - Complex *pwork = work.fortran_vec (); - int lrwork = (5 * (m < n ? m : n)) - 4; lrwork = lrwork > 1 ? lrwork : 1; Array<double> rwork (lrwork); double *prwork = rwork.fortran_vec (); + // Ask ZGELSS what the dimension of WORK should be. + + int lwork = -1; + + Array<Complex> work (1); + F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult, - nrr, ps, rcond, rank, pwork, lwork, - prwork, info)); + nrr, ps, rcond, rank, + work.fortran_vec (), lwork, prwork, + info)); if (f77_exception_encountered) (*current_liboctave_error_handler) ("unrecoverable error in zgelss"); else { - retval.resize (n, nrhs); - for (int j = 0; j < nrhs; j++) - for (int i = 0; i < n; i++) - retval.elem (i, j) = result.elem (i, j); + lwork = static_cast<int> (real (work(0))); + work.resize (lwork); + + F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult, + nrr, ps, rcond, rank, + work.fortran_vec (), lwork, + prwork, info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in zgelss"); + else + { + retval.resize (n, nrhs); + for (int j = 0; j < nrhs; j++) + for (int i = 0; i < n; i++) + retval.elem (i, j) = result.elem (i, j); + } } } @@ -1634,33 +1644,43 @@ double rcond = -1.0; - int lwork; - if (m < n) - lwork = 2*m + (nrhs > n ? nrhs : n); - else - lwork = 2*n + (nrhs > m ? nrhs : m); - - lwork *= 16; - - Array<Complex> work (lwork); - Complex *pwork = work.fortran_vec (); - int lrwork = (5 * (m < n ? m : n)) - 4; lrwork = lrwork > 1 ? lrwork : 1; Array<double> rwork (lrwork); double *prwork = rwork.fortran_vec (); + // Ask ZGELSS what the dimension of WORK should be. + + int lwork = -1; + + Array<Complex> work (1); + F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult, - nrr, ps, rcond, rank, pwork, lwork, - prwork, info)); + nrr, ps, rcond, rank, + work.fortran_vec (), lwork, prwork, + info)); if (f77_exception_encountered) (*current_liboctave_error_handler) ("unrecoverable error in zgelss"); else { - retval.resize (n); - for (int i = 0; i < n; i++) - retval.elem (i) = result.elem (i); + lwork = static_cast<int> (real (work(0))); + work.resize (lwork); + + F77_XFCN (zgelss, ZGELSS, (m, n, nrhs, tmp_data, m, presult, + nrr, ps, rcond, rank, + work.fortran_vec (), lwork, + prwork, info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in zgelss"); + else + { + retval.resize (n); + for (int i = 0; i < n; i++) + retval.elem (i) = result.elem (i); + } } }
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,8 @@ +2000-12-13 John W. Eaton <jwe@bevo.che.wisc.edu> + + * dMatrix.cc (Matrix::lssolve): Ask DGELSS for size of work vector. + * CMatrix.cc (ComplexMatrix::lssolve): Likewise, for ZGELSS. + 2000-12-09 John W. Eaton <jwe@bevo.che.wisc.edu> * Range.cc (Range::nelem_internal): Call round here, not tfloor.
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -1200,32 +1200,37 @@ double rcond = -1.0; - int lwork; - if (m < n) - lwork = 3*m + (2*m > nrhs - ? (2*m > n ? 2*m : n) - : (nrhs > n ? nrhs : n)); - else - lwork = 3*n + (2*n > nrhs - ? (2*n > m ? 2*n : m) - : (nrhs > m ? nrhs : m)); - - lwork *= 16; - - Array<double> work (lwork); - double *pwork = work.fortran_vec (); + // Ask DGELSS what the dimension of WORK should be. + + int lwork = -1; + + Array<double> work (1); F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult, nrr, ps, - rcond, rank, pwork, lwork, info)); + rcond, rank, work.fortran_vec (), + lwork, info)); if (f77_exception_encountered) (*current_liboctave_error_handler) ("unrecoverable error in dgelss"); else { - retval.resize (n, nrhs); - for (int j = 0; j < nrhs; j++) - for (int i = 0; i < n; i++) - retval.elem (i, j) = result.elem (i, j); + lwork = static_cast<int> (work(0)); + work.resize (lwork); + + F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult, + nrr, ps, rcond, rank, + work.fortran_vec (), lwork, info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dgelss"); + else + { + retval.resize (n, nrhs); + for (int j = 0; j < nrhs; j++) + for (int i = 0; i < n; i++) + retval.elem (i, j) = result.elem (i, j); + } } } @@ -1303,31 +1308,36 @@ double rcond = -1.0; - int lwork; - if (m < n) - lwork = 3*m + (2*m > nrhs - ? (2*m > n ? 2*m : n) - : (nrhs > n ? nrhs : n)); - else - lwork = 3*n + (2*n > nrhs - ? (2*n > m ? 2*n : m) - : (nrhs > m ? nrhs : m)); - - lwork *= 16; - - Array<double> work (lwork); - double *pwork = work.fortran_vec (); - - F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult, nrr, - ps, rcond, rank, pwork, lwork, info)); + // Ask DGELSS what the dimension of WORK should be. + + int lwork = -1; + + Array<double> work (1); + + F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult, nrr, ps, + rcond, rank, work.fortran_vec (), + lwork, info)); if (f77_exception_encountered) (*current_liboctave_error_handler) ("unrecoverable error in dgelss"); else { - retval.resize (n); - for (int i = 0; i < n; i++) - retval.elem (i) = result.elem (i); + lwork = static_cast<int> (work(0)); + work.resize (lwork); + + F77_XFCN (dgelss, DGELSS, (m, n, nrhs, tmp_data, m, presult, + nrr, ps, rcond, rank, + work.fortran_vec (), lwork, info)); + + if (f77_exception_encountered) + (*current_liboctave_error_handler) + ("unrecoverable error in dgelss"); + else + { + retval.resize (n); + for (int i = 0; i < n; i++) + retval.elem (i) = result.elem (i); + } } }