changeset 7544:f9983d2761df

more xGELSD workspace fixes
author Jaroslav Hajek <highegg@gmail.com>
date Mon, 03 Mar 2008 02:18:12 -0500
parents b84c5cbc0812
children 5b806195190d
files liboctave/CMatrix.cc liboctave/ChangeLog liboctave/dMatrix.cc
diffstat 3 files changed, 29 insertions(+), 14 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/CMatrix.cc	Fri Feb 29 04:09:03 2008 -0500
+++ b/liboctave/CMatrix.cc	Mon Mar 03 02:18:12 2008 -0500
@@ -2439,11 +2439,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;
 
@@ -2620,8 +2620,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
@@ -2629,11 +2633,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	Fri Feb 29 04:09:03 2008 -0500
+++ b/liboctave/ChangeLog	Mon Mar 03 02:18:12 2008 -0500
@@ -1,3 +1,10 @@
+2008-01-25  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-27  John W. Eaton  <jwe@octave.org>
 
 	* oct-rand.cc (class octave_rand): Make it a proper singleton class.
--- a/liboctave/dMatrix.cc	Fri Feb 29 04:09:03 2008 -0500
+++ b/liboctave/dMatrix.cc	Mon Mar 03 02:18:12 2008 -0500
@@ -2053,7 +2053,7 @@
 #else
       double tmp = log (dminmn / dsmlsizp1) / log (2.0);
 #endif
-      octave_idx_type nlvl = static_cast<int> (tmp) + 1;
+      octave_idx_type nlvl = static_cast<octave_idx_type> (tmp) + 1;
       if (nlvl < 0)
 	nlvl = 0;
 
@@ -2231,19 +2231,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;