Mercurial > octave-nkf

changeset 4725:fa612b2cbfe9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-01-23 16:42:51 by jwe]
author jwe
date Fri, 23 Jan 2004 16:42:51 +0000
parents bdacd0383fbd
children 14dc2267c343
files liboctave/ChangeLog liboctave/EIG.cc liboctave/EIG.h src/ChangeLog src/DLD-FUNCTIONS/eig.cc
diffstat 5 files changed, 46 insertions(+), 27 deletions(-) [+]
line wrap: on
line diff
--- a/liboctave/ChangeLog	Fri Jan 23 14:58:29 2004 +0000
+++ b/liboctave/ChangeLog	Fri Jan 23 16:42:51 2004 +0000
@@ -1,3 +1,10 @@
+2004-01-23  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* EIG.cc, EIG.h (EIG::init, EIG::symmetric_init, EIG::hermitian_init):
+	New arg, calc_eigenvectors.
+	* EIG.h (EIG:EIG): New optional arg, calc_eigenvectors.
+	Based on patch from David Bateman <dbateman@free.fr>.
+
 2004-01-22  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* Array.cc (Array<T>::assign2, Array<T>::assignN):
--- a/liboctave/EIG.cc	Fri Jan 23 14:58:29 2004 +0000
+++ b/liboctave/EIG.cc	Fri Jan 23 16:42:51 2004 +0000
@@ -71,10 +71,10 @@
 }
 
 int
-EIG::init (const Matrix& a)
+EIG::init (const Matrix& a, bool calc_ev)
 {
   if (a.is_symmetric ())
-    return symmetric_init (a);
+    return symmetric_init (a, calc_ev);
 
   int n = a.rows ();
 
@@ -95,7 +95,8 @@
   Array<double> wi (n);
   double *pwi = wi.fortran_vec ();
 
-  Matrix vr (n, n);
+  int nvr = calc_ev ? n : 0;
+  Matrix vr (nvr, nvr);
   double *pvr = vr.fortran_vec ();
 
   // XXX FIXME XXX -- it might be possible to choose a better value of
@@ -109,7 +110,7 @@
   int idummy = 1;
 
   F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-			   F77_CONST_CHAR_ARG2 ("V", 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)
@@ -124,14 +125,14 @@
       else
 	{
 	  lambda.resize (n);
-	  v.resize (n, n);
+	  v.resize (nvr, nvr);
 
 	  for (int j = 0; j < n; j++)
 	    {
 	      if (wi.elem (j) == 0.0)
 		{
 		  lambda.elem (j) = Complex (wr.elem (j));
-		  for (int i = 0; i < n; i++)
+		  for (int i = 0; i < nvr; i++)
 		    v.elem (i, j) = vr.elem (i, j);
 		}
 	      else
@@ -146,7 +147,7 @@
 		  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 < n; i++)
+		  for (int i = 0; i < nvr; i++)
 		    {
 		      double real_part = vr.elem (i, j);
 		      double imag_part = vr.elem (i, j+1);
@@ -163,7 +164,7 @@
 }
 
 int
-EIG::symmetric_init (const Matrix& a)
+EIG::symmetric_init (const Matrix& a, bool calc_ev)
 {
   int n = a.rows ();
 
@@ -188,7 +189,7 @@
   Array<double> work (lwork);
   double *pwork = work.fortran_vec ();
 
-  F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 ("V", 1),
+  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)
@@ -201,17 +202,17 @@
   else
     {
       lambda = ComplexColumnVector (wr);
-      v = ComplexMatrix (atmp);
+      v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
     }
 
   return info;
 }
 
 int
-EIG::init (const ComplexMatrix& a)
+EIG::init (const ComplexMatrix& a, bool calc_ev)
 {
   if (a.is_hermitian ())
-    return hermitian_init (a);
+    return hermitian_init (a, calc_ev);
 
   int n = a.rows ();
 
@@ -229,7 +230,8 @@
   ComplexColumnVector w (n);
   Complex *pw = w.fortran_vec ();
 
-  ComplexMatrix vtmp (n, n);
+  int nvr = calc_ev ? n : 0;
+  ComplexMatrix vtmp (nvr, nvr);
   Complex *pv = vtmp.fortran_vec ();
 
   // XXX FIXME XXX -- it might be possible to choose a better value of
@@ -247,7 +249,7 @@
   int idummy = 1;
 
   F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
-			   F77_CONST_CHAR_ARG2 ("V", 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)
@@ -267,7 +269,7 @@
 }
 
 int
-EIG::hermitian_init (const ComplexMatrix& a)
+EIG::hermitian_init (const ComplexMatrix& a, bool calc_ev)
 {
   int n = a.rows ();
 
@@ -296,7 +298,7 @@
   Array<double> rwork (lrwork);
   double *prwork = rwork.fortran_vec ();
 
-  F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 ("V", 1),
+  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)
@@ -309,7 +311,7 @@
   else
     {
       lambda = ComplexColumnVector (wr);
-      v = ComplexMatrix (atmp);
+      v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
     }
 
   return info;
--- a/liboctave/EIG.h	Fri Jan 23 14:58:29 2004 +0000
+++ b/liboctave/EIG.h	Fri Jan 23 16:42:51 2004 +0000
@@ -44,13 +44,17 @@
   EIG (void)
     : lambda (), v () { }
 
-  EIG (const Matrix& a) { init (a); }
+  EIG (const Matrix& a, bool calc_eigenvectors = true)
+    { init (a, calc_eigenvectors); }
 
-  EIG (const Matrix& a, int& info) { info = init (a); }
+  EIG (const Matrix& a, int& info, bool calc_eigenvectors = true)
+    { info = init (a, calc_eigenvectors); }
 
-  EIG (const ComplexMatrix& a) { init (a); }
+  EIG (const ComplexMatrix& a, bool calc_eigenvectors = true)
+    { init (a, calc_eigenvectors); }
 
-  EIG (const ComplexMatrix& a, int& info) { info = init (a); }
+  EIG (const ComplexMatrix& a, int& info, bool calc_eigenvectors = true)
+    { info = init (a, calc_eigenvectors); }
 
   EIG (const EIG& a)
     : lambda (a.lambda), v (a.v) { }
@@ -78,11 +82,11 @@
   ComplexColumnVector lambda;
   ComplexMatrix v;
 
-  int init (const Matrix& a);
-  int init (const ComplexMatrix& a);
+  int init (const Matrix& a, bool calc_eigenvectors);
+  int init (const ComplexMatrix& a, bool calc_eigenvectors);
 
-  int symmetric_init (const Matrix& a);
-  int hermitian_init (const ComplexMatrix& a);
+  int symmetric_init (const Matrix& a, bool calc_eigenvectors);
+  int hermitian_init (const ComplexMatrix& a, bool calc_eigenvectors);
 };
 
 #endif
--- a/src/ChangeLog	Fri Jan 23 14:58:29 2004 +0000
+++ b/src/ChangeLog	Fri Jan 23 16:42:51 2004 +0000
@@ -1,3 +1,9 @@
+2004-01-23  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* DLD-FUNCTIONS/eig.cc (Feig): Use new optional arg for EIG to
+	avoid computing eigenvectors if not requested.
+	Based on a patch from David Bateman  <dbateman@free.fr>.
+
 2004-01-23  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* ov-cell.cc (all_strings): Always compute total required length
--- a/src/DLD-FUNCTIONS/eig.cc	Fri Jan 23 14:58:29 2004 +0000
+++ b/src/DLD-FUNCTIONS/eig.cc	Fri Jan 23 16:42:51 2004 +0000
@@ -81,7 +81,7 @@
       if (error_state)
 	return retval;
       else
-	result = EIG (tmp);
+	result = EIG (tmp, nargout > 1);
     }
   else if (arg.is_complex_type ())
     {
@@ -90,7 +90,7 @@
       if (error_state)
 	return retval;
       else
-	result = EIG (ctmp);
+	result = EIG (ctmp, nargout > 1);
     }
   else
     {