Mercurial > hg > octave-lyh
changeset 7071:c3b479e753dd
[project @ 2007-10-26 15:14:34 by jwe]
| author | jwe |
|---|---|
| date | Fri, 26 Oct 2007 15:14:35 +0000 |
| parents | 7593f8e83a2e |
| children | b48d486f641d |
| files | PROJECTS doc/interpreter/linalg.txi liboctave/CMatrix.cc liboctave/ChangeLog liboctave/dMatrix.cc |
| diffstat | 5 files changed, 55 insertions(+), 42 deletions(-) [+] |
line wrap: on
line diff
--- a/PROJECTS +++ b/PROJECTS @@ -66,7 +66,7 @@ * Consider making the behavior of the / and \ operators for non-square systems compatible with Matlab. Currently, they return - the minimum norm solution from DGELSY, which behaves differently. + the minimum norm solution from DGELSS, which behaves differently. --------------- Sparse Matrices:
--- a/doc/interpreter/linalg.txi +++ b/doc/interpreter/linalg.txi @@ -63,7 +63,7 @@ @item If the matrix is not square, or any of the previous solvers flags a singular or near singular matrix, find a least squares solution using -the @sc{Lapack} xGELSY function. +the @sc{Lapack} xGELSS function. @end enumerate The user can force the type of the matrix with the @code{matrix_type}
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -120,9 +120,9 @@ F77_CHAR_ARG_LEN_DECL); F77_RET_T - F77_FUNC (zgelsy, ZGELSY) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + F77_FUNC (zgelss, ZGELSS) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, Complex*, const octave_idx_type&, Complex*, - const octave_idx_type&, octave_idx_type*, double&, octave_idx_type&, + const octave_idx_type&, double*, double&, octave_idx_type&, Complex*, const octave_idx_type&, double*, octave_idx_type&); F77_RET_T @@ -2448,40 +2448,43 @@ Complex *presult = result.fortran_vec (); - Array<octave_idx_type> jpvt (n); - octave_idx_type *pjpvt = jpvt.fortran_vec (); + octave_idx_type len_s = m < n ? m : n; + Array<double> s (len_s); + double *ps = s.fortran_vec (); double rcond = -1.0; - Array<double> rwork (2 * n); + octave_idx_type lrwork = (5 * (m < n ? m : n)) - 4; + lrwork = lrwork > 1 ? lrwork : 1; + Array<double> rwork (lrwork); double *prwork = rwork.fortran_vec (); - // Ask ZGELSY what the dimension of WORK should be. + // Ask ZGELSS what the dimension of WORK should be. octave_idx_type lwork = -1; Array<Complex> work (1); - F77_XFCN (zgelsy, ZGELSY, (m, n, nrhs, tmp_data, m, presult, - nrr, pjpvt, rcond, rank, + 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 zgelsy"); + (*current_liboctave_error_handler) ("unrecoverable error in zgelss"); else { lwork = static_cast<octave_idx_type> (std::real (work(0))); work.resize (lwork); - F77_XFCN (zgelsy, ZGELSY, (m, n, nrhs, tmp_data, m, presult, - nrr, pjpvt, rcond, rank, + 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 zgelsy"); + ("unrecoverable error in zgelss"); else { retval.resize (n, nrhs); @@ -2560,40 +2563,43 @@ Complex *presult = result.fortran_vec (); - Array<octave_idx_type> jpvt (n); - octave_idx_type *pjpvt = jpvt.fortran_vec (); + octave_idx_type len_s = m < n ? m : n; + Array<double> s (len_s); + double *ps = s.fortran_vec (); double rcond = -1.0; - Array<double> rwork (2 * n); + octave_idx_type lrwork = (5 * (m < n ? m : n)) - 4; + lrwork = lrwork > 1 ? lrwork : 1; + Array<double> rwork (lrwork); double *prwork = rwork.fortran_vec (); - // Ask ZGELSY what the dimension of WORK should be. + // Ask ZGELSS what the dimension of WORK should be. octave_idx_type lwork = -1; Array<Complex> work (1); - F77_XFCN (zgelsy, ZGELSY, (m, n, nrhs, tmp_data, m, presult, - nrr, pjpvt, rcond, rank, + 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 zgelsy"); + (*current_liboctave_error_handler) ("unrecoverable error in zgelss"); else { lwork = static_cast<int> (std::real (work(0))); work.resize (lwork); - F77_XFCN (zgelsy, ZGELSY, (m, n, nrhs, tmp_data, m, presult, - nrr, pjpvt, rcond, rank, + 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 zgelsy"); + ("unrecoverable error in zgelss"); else { retval.resize (n);
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,8 @@ +2007-09-26 John W. Eaton <jwe@octave.org> + + * dMatrix.cc (lssolve): Revert change of 2007-09-26. + * CMatrix.cc (lssolve): Ditto. + 2007-10-25 John W. Eaton <jwe@octave.org> * oct-time.cc (octave_gmtime::init, octave_localtime::init):
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -117,9 +117,9 @@ F77_CHAR_ARG_LEN_DECL); F77_RET_T - F77_FUNC (dgelsy, DGELSY) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, + F77_FUNC (dgelss, DGELSS) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, double*, const octave_idx_type&, double*, - const octave_idx_type&, octave_idx_type*, double&, octave_idx_type&, + const octave_idx_type&, double*, double&, octave_idx_type&, double*, const octave_idx_type&, octave_idx_type&); F77_RET_T @@ -2072,35 +2072,36 @@ double *presult = result.fortran_vec (); - Array<octave_idx_type> jpvt (n); - octave_idx_type *pjpvt = jpvt.fortran_vec (); + octave_idx_type len_s = m < n ? m : n; + Array<double> s (len_s); + double *ps = s.fortran_vec (); double rcond = -1.0; - // Ask DGELSY what the dimension of WORK should be. + // Ask DGELSS what the dimension of WORK should be. octave_idx_type lwork = -1; Array<double> work (1); - F77_XFCN (dgelsy, DGELSY, (m, n, nrhs, tmp_data, m, presult, nrr, pjpvt, + 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 dgelsy"); + (*current_liboctave_error_handler) ("unrecoverable error in dgelss"); else { lwork = static_cast<octave_idx_type> (work(0)); work.resize (lwork); - F77_XFCN (dgelsy, DGELSY, (m, n, nrhs, tmp_data, m, presult, - nrr, pjpvt, rcond, rank, + 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 dgelsy"); + ("unrecoverable error in dgelss"); else { retval.resize (n, nrhs); @@ -2181,35 +2182,36 @@ double *presult = result.fortran_vec (); - Array<octave_idx_type> jpvt (n); - octave_idx_type *pjpvt = jpvt.fortran_vec (); + octave_idx_type len_s = m < n ? m : n; + Array<double> s (len_s); + double *ps = s.fortran_vec (); double rcond = -1.0; - // Ask DGELSY what the dimension of WORK should be. + // Ask DGELSS what the dimension of WORK should be. octave_idx_type lwork = -1; Array<double> work (1); - F77_XFCN (dgelsy, DGELSY, (m, n, nrhs, tmp_data, m, presult, nrr, pjpvt, + 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 dgelsy"); + (*current_liboctave_error_handler) ("unrecoverable error in dgelss"); else { lwork = static_cast<octave_idx_type> (work(0)); work.resize (lwork); - F77_XFCN (dgelsy, DGELSY, (m, n, nrhs, tmp_data, m, presult, - nrr, pjpvt, rcond, rank, + 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 dgelsy"); + ("unrecoverable error in dgelss"); else { retval.resize (n);