Mercurial > octave-nkf

changeset 7071:c3b479e753dd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[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	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:
--- 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}
--- 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<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	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  <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	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<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);