Mercurial > octave-nkf

changeset 4800:c322edde72ac

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-02-25 03:54:33 by jwe]
author jwe
date Wed, 25 Feb 2004 03:54:33 +0000
parents e2d7d1ef5e55
children b022780ac0b4
files liboctave/ChangeLog liboctave/EIG.cc
diffstat 2 files changed, 142 insertions(+), 74 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/ChangeLog	Tue Feb 24 19:47:18 2004 +0000
+++ b/liboctave/ChangeLog	Wed Feb 25 03:54:33 2004 +0000
@@ -1,3 +1,8 @@
+2004-02-24  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* EIG.cc (EIG::init, EIG::symmetric_init):
+	Query Lapack for workspace size.
+
 2004-02-23  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* Array.cc (Array<T>::resize_and_fill (const dim_vector&, const T&)): 
--- a/liboctave/EIG.cc	Tue Feb 24 19:47:18 2004 +0000
+++ b/liboctave/EIG.cc	Wed Feb 25 03:54:33 2004 +0000
@@ -99,12 +99,8 @@
   Matrix vr (nvr, nvr);
   double *pvr = vr.fortran_vec ();
 
-  // XXX FIXME XXX -- it might be possible to choose a better value of
-  // lwork that would result in more efficient computations.
-
-  int lwork = 8*n;
-  Array<double> work (lwork);
-  double *pwork = work.fortran_vec ();
+  int lwork = -1;
+  double dummy_work;
 
   double *dummy = 0;
   int idummy = 1;
@@ -112,53 +108,70 @@
   F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
 			   F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
 			   n, tmp_data, n, pwr, pwi, dummy,
-			   idummy, pvr, n, pwork, lwork, info
+			   idummy, pvr, n, &dummy_work, lwork, info
 			   F77_CHAR_ARG_LEN (1)
 			   F77_CHAR_ARG_LEN (1)));
 
-  if (f77_exception_encountered || info < 0)
-    (*current_liboctave_error_handler) ("unrecoverable error in dgeev");
-  else
+  if (! f77_exception_encountered && info == 0)
     {
+      lwork = static_cast<int> (dummy_work);
+      Array<double> work (lwork);
+      double *pwork = work.fortran_vec ();
+
+      F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			       F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       n, tmp_data, n, pwr, pwi, dummy,
+			       idummy, pvr, n, pwork, lwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+
+      if (f77_exception_encountered || info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in dgeev");
+	  return info;
+	}
+
       if (info > 0)
-	(*current_liboctave_error_handler) ("dgeev failed to converge");
-      else
 	{
-	  lambda.resize (n);
-	  v.resize (nvr, nvr);
+	  (*current_liboctave_error_handler) ("dgeev failed to converge");
+	  return info;
+	}
 
-	  for (int j = 0; j < n; j++)
+      lambda.resize (n);
+      v.resize (nvr, nvr);
+
+      for (int j = 0; j < n; j++)
+	{
+	  if (wi.elem (j) == 0.0)
 	    {
-	      if (wi.elem (j) == 0.0)
+	      lambda.elem (j) = Complex (wr.elem (j));
+	      for (int i = 0; i < nvr; i++)
+		v.elem (i, j) = vr.elem (i, j);
+	    }
+	  else
+	    {
+	      if (j+1 >= n)
 		{
-		  lambda.elem (j) = Complex (wr.elem (j));
-		  for (int i = 0; i < nvr; i++)
-		    v.elem (i, j) = vr.elem (i, j);
+		  (*current_liboctave_error_handler) ("EIG: internal error");
+		  return -1;
 		}
-	      else
-		{
-		  if (j+1 >= n)
-		    {
-		      (*current_liboctave_error_handler)
-			("EIG: internal error");
-		      return -1;
-		    }
+
+	      lambda.elem(j) = Complex (wr.elem(j), wi.elem(j));
+	      lambda.elem(j+1) = Complex (wr.elem(j+1), wi.elem(j+1));
 
-		  lambda.elem(j) = Complex (wr.elem(j), wi.elem(j));
-		  lambda.elem(j+1) = Complex (wr.elem(j+1), wi.elem(j+1));
-
-		  for (int i = 0; i < nvr; i++)
-		    {
-		      double real_part = vr.elem (i, j);
-		      double imag_part = vr.elem (i, j+1);
-		      v.elem (i, j) = Complex (real_part, imag_part);
-		      v.elem (i, j+1) = Complex (real_part, -imag_part);
-		    }
-		  j++;
+	      for (int i = 0; i < nvr; i++)
+		{
+		  double real_part = vr.elem (i, j);
+		  double imag_part = vr.elem (i, j+1);
+		  v.elem (i, j) = Complex (real_part, imag_part);
+		  v.elem (i, j+1) = Complex (real_part, -imag_part);
 		}
+	      j++;
 	    }
 	}
     }
+  else
+    (*current_liboctave_error_handler) ("dgeev workspace query failed");
 
   return info;
 }
@@ -182,28 +195,44 @@
   ColumnVector wr (n);
   double *pwr = wr.fortran_vec ();
 
-  // XXX FIXME XXX -- it might be possible to choose a better value of
-  // lwork that would result in more efficient computations.
-
-  int lwork = 8*n;
-  Array<double> work (lwork);
-  double *pwork = work.fortran_vec ();
+  int lwork = -1;
+  double dummy_work;
 
   F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
 			   F77_CONST_CHAR_ARG2 ("U", 1),
-			   n, tmp_data, n, pwr, pwork, lwork, info
+			   n, tmp_data, n, pwr, &dummy_work, lwork, info
 			   F77_CHAR_ARG_LEN (1)
 			   F77_CHAR_ARG_LEN (1)));
 
-  if (f77_exception_encountered || info < 0)
-    (*current_liboctave_error_handler) ("unrecoverable error in dsyev");
-  else if (info > 0)
-    (*current_liboctave_error_handler) ("dsyev failed to converge");
-  else
+  if (! f77_exception_encountered && info == 0)
     {
+      lwork = static_cast<int> (dummy_work);
+      Array<double> work (lwork);
+      double *pwork = work.fortran_vec ();
+
+      F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       F77_CONST_CHAR_ARG2 ("U", 1),
+			       n, tmp_data, n, pwr, pwork, lwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+
+      if (f77_exception_encountered || info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in dsyev");
+	  return info;
+	}
+
+      if (info > 0)
+	{
+	  (*current_liboctave_error_handler) ("dsyev failed to converge");
+	  return info;
+	}
+
       lambda = ComplexColumnVector (wr);
       v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
     }
+  else
+    (*current_liboctave_error_handler) ("dsyev workspace query failed");
 
   return info;
 }
@@ -234,12 +263,8 @@
   ComplexMatrix vtmp (nvr, nvr);
   Complex *pv = vtmp.fortran_vec ();
 
-  // XXX FIXME XXX -- it might be possible to choose a better value of
-  // lwork that would result in more efficient computations.
-
-  int lwork = 8*n;
-  Array<Complex> work (lwork);
-  Complex *pwork = work.fortran_vec ();
+  int lwork = -1;
+  Complex dummy_work;
 
   int lrwork = 2*n;
   Array<double> rwork (lrwork);
@@ -251,19 +276,40 @@
   F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
 			   F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
 			   n, tmp_data, n, pw, dummy, idummy,
-			   pv, n, pwork, lwork, prwork, info
+			   pv, n, &dummy_work, lwork, prwork, info
 			   F77_CHAR_ARG_LEN (1)
 			   F77_CHAR_ARG_LEN (1)));
 
-  if (f77_exception_encountered || info < 0)
-    (*current_liboctave_error_handler) ("unrecoverable error in zgeev");
-  else if (info > 0)
-    (*current_liboctave_error_handler) ("zgeev failed to converge");
-  else
+  if (! f77_exception_encountered && info == 0)
     {
+      lwork = static_cast<int> (dummy_work.real ());
+      Array<Complex> work (lwork);
+      Complex *pwork = work.fortran_vec ();
+
+      F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			       F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       n, tmp_data, n, pw, dummy, idummy,
+			       pv, n, pwork, lwork, prwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+
+      if (f77_exception_encountered || info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in zgeev");
+	  return info;
+	}
+
+      if (info > 0)
+	{
+	  (*current_liboctave_error_handler) ("zgeev failed to converge");
+	  return info;
+	}
+
       lambda = w;
       v = vtmp;
     }
+  else
+    (*current_liboctave_error_handler) ("zgeev workspace query failed");
 
   return info;
 }
@@ -287,12 +333,8 @@
   ColumnVector wr (n);
   double *pwr = wr.fortran_vec ();
 
-  // XXX FIXME XXX -- it might be possible to choose a better value of
-  // lwork that would result in more efficient computations.
-
-  int lwork = 8*n;
-  Array<Complex> work (lwork);
-  Complex *pwork = work.fortran_vec ();
+  int lwork = -1;
+  Complex dummy_work;
 
   int lrwork = 3*n;
   Array<double> rwork (lrwork);
@@ -300,19 +342,40 @@
 
   F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
 			   F77_CONST_CHAR_ARG2 ("U", 1),
-			   n, tmp_data, n, pwr, pwork, lwork, prwork, info
+			   n, tmp_data, n, pwr, &dummy_work, lwork,
+			   prwork, info
 			   F77_CHAR_ARG_LEN (1)
 			   F77_CHAR_ARG_LEN (1)));
 
-  if (f77_exception_encountered || info < 0)
-    (*current_liboctave_error_handler) ("unrecoverable error in zheev");
-  else if (info > 0)
-    (*current_liboctave_error_handler) ("zheev failed to converge");
-  else
+  if (! f77_exception_encountered && info == 0)
     {
+      lwork = static_cast<int> (dummy_work.real ());
+      Array<Complex> work (lwork);
+      Complex *pwork = work.fortran_vec ();
+
+      F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       F77_CONST_CHAR_ARG2 ("U", 1),
+			       n, tmp_data, n, pwr, pwork, lwork, prwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+
+      if (f77_exception_encountered || info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in zheev");
+	  return info;
+	}
+
+      if (info > 0)
+	{
+	  (*current_liboctave_error_handler) ("zheev failed to converge");
+	  return info;
+	}
+
       lambda = ComplexColumnVector (wr);
       v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
     }
+  else
+    (*current_liboctave_error_handler) ("zheev workspace query failed");
 
   return info;
 }