# HG changeset patch # User jwe # Date 1193411675 0 # Node ID c3b479e753dd21d1f8585aa530be5f9077a7c181 # Parent 7593f8e83a2e52704f886fac38262ca834df9f05 [project @ 2007-10-26 15:14:34 by jwe] diff -r 7593f8e83a2e -r c3b479e753dd PROJECTS --- a/PROJECTS Thu Oct 25 20:41:17 2007 +0000 +++ b/PROJECTS Fri Oct 26 15:14:35 2007 +0000 @@ -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: diff -r 7593f8e83a2e -r c3b479e753dd doc/interpreter/linalg.txi --- a/doc/interpreter/linalg.txi Thu Oct 25 20:41:17 2007 +0000 +++ b/doc/interpreter/linalg.txi Fri Oct 26 15:14:35 2007 +0000 @@ -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} diff -r 7593f8e83a2e -r c3b479e753dd liboctave/CMatrix.cc --- a/liboctave/CMatrix.cc Thu Oct 25 20:41:17 2007 +0000 +++ b/liboctave/CMatrix.cc Fri Oct 26 15:14:35 2007 +0000 @@ -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 jpvt (n); - octave_idx_type *pjpvt = jpvt.fortran_vec (); + octave_idx_type len_s = m < n ? m : n; + Array s (len_s); + double *ps = s.fortran_vec (); double rcond = -1.0; - Array rwork (2 * n); + octave_idx_type lrwork = (5 * (m < n ? m : n)) - 4; + lrwork = lrwork > 1 ? lrwork : 1; + Array 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 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 (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 jpvt (n); - octave_idx_type *pjpvt = jpvt.fortran_vec (); + octave_idx_type len_s = m < n ? m : n; + Array s (len_s); + double *ps = s.fortran_vec (); double rcond = -1.0; - Array rwork (2 * n); + octave_idx_type lrwork = (5 * (m < n ? m : n)) - 4; + lrwork = lrwork > 1 ? lrwork : 1; + Array 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 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 (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); diff -r 7593f8e83a2e -r c3b479e753dd liboctave/ChangeLog --- a/liboctave/ChangeLog Thu Oct 25 20:41:17 2007 +0000 +++ b/liboctave/ChangeLog Fri Oct 26 15:14:35 2007 +0000 @@ -1,3 +1,8 @@ +2007-09-26 John W. Eaton + + * dMatrix.cc (lssolve): Revert change of 2007-09-26. + * CMatrix.cc (lssolve): Ditto. + 2007-10-25 John W. Eaton * oct-time.cc (octave_gmtime::init, octave_localtime::init): diff -r 7593f8e83a2e -r c3b479e753dd liboctave/dMatrix.cc --- a/liboctave/dMatrix.cc Thu Oct 25 20:41:17 2007 +0000 +++ b/liboctave/dMatrix.cc Fri Oct 26 15:14:35 2007 +0000 @@ -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 jpvt (n); - octave_idx_type *pjpvt = jpvt.fortran_vec (); + octave_idx_type len_s = m < n ? m : n; + Array 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 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 (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 jpvt (n); - octave_idx_type *pjpvt = jpvt.fortran_vec (); + octave_idx_type len_s = m < n ? m : n; + Array 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 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 (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);