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	Sun Dec 10 06:03:06 2000 +0000
+++ b/liboctave/CMatrix.cc	Wed Dec 13 19:02:43 2000 +0000
@@ -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	Sun Dec 10 06:03:06 2000 +0000
+++ b/liboctave/ChangeLog	Wed Dec 13 19:02:43 2000 +0000
@@ -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	Sun Dec 10 06:03:06 2000 +0000
+++ b/liboctave/dMatrix.cc	Wed Dec 13 19:02:43 2000 +0000
@@ -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);
+	    }
 	}
     }