From octave-maintainers-request at bevo dot che dot wisc dot edu Fri Jan 23 10:44:17 2004 Subject: eigenvalues 3 times speedup patch [Was: benchmarks] From: "John W. Eaton" To: David Bateman Cc: Paul Kienzle , octave-maintainers@bevo.che.wisc.edu Date: Fri, 23 Jan 2004 10:42:28 -0600 On 23-Jan-2004, David Bateman 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 at @ -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 (); at @ -95,7 +95,8 @@ Array 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 at @ -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) at @ -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 at @ -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); at @ -163,7 +164,7 @@ } int -EIG::symmetric_init (const Matrix& a) +EIG::symmetric_init (const Matrix& a, bool calc_ev) { int n = a.rows (); at @ -188,7 +189,7 @@ Array 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) at @ -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 (); at @ -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 at @ -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) at @ -267,7 +269,7 @@ } int -EIG::hermitian_init (const ComplexMatrix& a) +EIG::hermitian_init (const ComplexMatrix& a, bool calc_ev) { int n = a.rows (); at @ -296,7 +298,7 @@ Array 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) at @ -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 at @ -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) { } at @ -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 at @ -81,7 +81,7 @@ if (error_state) return retval; else - result = EIG (tmp); + result = EIG (tmp, nargout > 1); } else if (arg.is_complex_type ()) { at @ -90,7 +90,7 @@ if (error_state) return retval; else - result = EIG (ctmp); + result = EIG (ctmp, nargout > 1); } else {