eigenvalues 3 times speedup patch [Was: benchmarks]

[John W. Eaton]

On 23-Jan-2004, David Bateman <address@hidden> wrote:

| Ok, this e-mail of your has drawn a lot of replies from me. I think this
| is probably the last thread I'll draw from it for now. However, check 
| 
| > II.B Matlab 0.86 - Octave 2.30 - O-Matrix 0.44
| >     b = eig(a);
| 
| Well, there is a simple reason why octave is slower here, for this benchmark.
| The benchmark asks for the eigenvalues only, but octave calculates the
| eignevectors as well... Dohhhh....
| 
| I've written a small patch to 2.1.53, that makes "eig(a)" only calculate the
| eigenvectors when asked for them.

OK, but how about the following patch instead?

It might also be good to fix these functions to use the optimal
workspace returned from Lapack.

Thanks,

jwe


Index: liboctave/EIG.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/EIG.cc,v
retrieving revision 1.24
diff -u -r1.24 EIG.cc
--- liboctave/EIG.cc    27 Oct 2003 20:38:02 -0000      1.24
+++ liboctave/EIG.cc    23 Jan 2004 16:26:07 -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;
Index: liboctave/EIG.h
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/EIG.h,v
retrieving revision 1.16
diff -u -r1.16 EIG.h
--- liboctave/EIG.h     20 Nov 2002 16:56:48 -0000      1.16
+++ liboctave/EIG.h     23 Jan 2004 16:26:07 -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
Index: src/DLD-FUNCTIONS/eig.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/DLD-FUNCTIONS/eig.cc,v
retrieving revision 1.5
diff -u -r1.5 eig.cc
--- src/DLD-FUNCTIONS/eig.cc    25 Nov 2003 05:41:36 -0000      1.5
+++ src/DLD-FUNCTIONS/eig.cc    23 Jan 2004 16:26:17 -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
     {