Mercurial > hg > octave-max

changeset 2815:33486d9e2d00

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1997-03-14 08:24:46 by jwe]
author jwe
date Fri, 14 Mar 1997 08:25:58 +0000
parents ffa60dc8e49b
children ad4bf2a82b4f
files NEWS kpathsea/ChangeLog kpathsea/configure.in libcruft/ChangeLog liboctave/Array2.cc liboctave/Array2.h liboctave/CMatrix.cc liboctave/CMatrix.h liboctave/ChangeLog liboctave/EIG.cc liboctave/EIG.h readline/ChangeLog readline/configure.in
diffstat 13 files changed, 250 insertions(+), 40 deletions(-) [+]
line wrap: on
line diff
--- a/NEWS
+++ b/NEWS
@@ -11,6 +11,10 @@
       function y = f (x) y = sqrt (x); endfunction
       quad ("f", 0, 1)
 
+  * If the argument to eig() is symmetric, Octave uses the specialized
+    Lapack subroutine for symmetric matrices for a significant
+    increase in performance.
+
 Summary of changes for version 2.0.5:
 ------------------------------------
 
--- a/kpathsea/ChangeLog
+++ b/kpathsea/ChangeLog
@@ -1,3 +1,8 @@
+Thu Mar 13 13:08:05 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* configure.in: Don't set special PIC options on SGI systems,
+	since PIC is apparently the default.
+
 Wed Mar 12 17:00:16 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* Makefile.in (install-strip): New target.
--- a/kpathsea/configure.in
+++ b/kpathsea/configure.in
@@ -37,6 +37,9 @@
     SHLEXT=sl
     SH_LDFLAGS="-shared -fPIC"
   ;;
+  *-sgi-*)
+    CPICFLAG=
+  ;;
 esac
 AC_MSG_RESULT([defining CPICFLAG to be $CPICFLAG])
 AC_MSG_RESULT([defining SHLEXT to be $SHLEXT])
--- a/libcruft/ChangeLog
+++ b/libcruft/ChangeLog
@@ -1,3 +1,8 @@
+Thu Mar 13 22:31:35 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* blas, lapack: Add new files for symmetric eigenvalue
+	computation.
+
 Wed Mar 12 16:59:59 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* misc/Makefile.in (install-strip): New target.
--- a/liboctave/Array2.cc
+++ b/liboctave/Array2.cc
@@ -171,6 +171,23 @@
   return *this;
 }
 
+template <class T>
+bool
+Array2<T>::is_symmetric (void) const
+{
+  if (is_square () && d1 > 0)
+    {
+      for (int i = 0; i < d1; i++)
+	for (int j = i+1; j < d2; j++)
+	  if (elem (i, j) != elem (j, i))
+	    return false;
+
+      return true;
+    }
+
+  return false;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
--- a/liboctave/Array2.h
+++ b/liboctave/Array2.h
@@ -173,6 +173,10 @@
 
   Array2<T>& insert (const Array2<T>& a, int r, int c);
 
+  bool is_square (void) const { return (d1 == d2); }
+
+  bool is_symmetric (void) const;
+
 #ifdef HEAVYWEIGHT_INDEXING
   void maybe_delete_elements (idx_vector& i, idx_vector& j);
 
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -181,6 +181,25 @@
   return !(*this == a);
 }
 
+bool
+ComplexMatrix::is_hermitian (void) const
+{
+  int nr = rows ();
+  int nc = cols ();
+
+  if (is_square () && nr > 0)
+    {
+      for (int i = 0; i < nr; i++)
+	for (int j = i; j < nc; j++)
+	  if (elem (i, j) != conj (elem (j, i)))
+	    return false;
+
+      return true;
+    }
+
+  return false;
+}
+
 // destructive insert/delete/reorder operations
 
 ComplexMatrix&
--- a/liboctave/CMatrix.h
+++ b/liboctave/CMatrix.h
@@ -72,6 +72,8 @@
   bool operator == (const ComplexMatrix& a) const;
   bool operator != (const ComplexMatrix& a) const;
 
+  bool is_hermitian (void) const;
+
   // destructive insert/delete/reorder operations
 
   ComplexMatrix& insert (const Matrix& a, int r, int c);
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,17 @@
+Fri Mar 14 00:29:46 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* EIG.cc (EIG::hermitian_init (const ComplexMatrix&)): New function.
+	(EIG::init (const ComplexMatrix&)): Call it if arg is hermitian.
+	(EIG::symmetric_init (const Matrix&)): New function.
+	(EIG::init (const Matrix&)): Call it if arg is symmetric.
+
+	* CMatrix.cc (ComplexMatrix::is_hermitian): New function.
+
+Thu Mar 13 17:04:26 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* Array2.cc (is_symmetric): New function.
+	* Array2.h (is_square): New function.
+
 Wed Mar 12 16:59:49 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* Makefile.in (install-strip): New target.
--- a/liboctave/EIG.cc
+++ b/liboctave/EIG.cc
@@ -29,6 +29,7 @@
 #endif
 
 #include "EIG.h"
+#include "dColVector.h"
 #include "f77-fcn.h"
 #include "lo-error.h"
 
@@ -45,11 +46,22 @@
 			      Complex*, const int&, Complex*,
 			      const int&, Complex*, const int&,
 			      double*, int&, long, long);
+
+  int F77_FCN (dsyev, DSYEV) (const char*, const char*, const int&,
+			      double*, const int&, double*, double*,
+			      const int&, int&, long, long);
+
+  int F77_FCN (zheev, ZHEEV) (const char*, const char*, const int&,
+			      Complex*, const int&, double*, Complex*,
+			      const int&, double*, int&, long, long);
 }
 
 int
 EIG::init (const Matrix& a)
 {
+  if (a.is_symmetric ())
+    return symmetric_init (a);
+
   int n = a.rows ();
 
   if (n != a.cols ())
@@ -58,7 +70,7 @@
       return -1;
     }
 
-  int info;
+  int info = 0;
 
   Matrix atmp = a;
   double *tmp_data = atmp.fortran_vec ();
@@ -72,8 +84,10 @@
   Matrix vr (n, n);
   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 ();
 
@@ -83,40 +97,95 @@
   F77_XFCN (dgeev, DGEEV, ("N", "V", n, tmp_data, n, pwr, pwi, dummy,
 			   idummy, pvr, n, pwork, lwork, info, 1L, 1L));
 
-  if (f77_exception_encountered)
+  if (f77_exception_encountered || info < 0)
     (*current_liboctave_error_handler) ("unrecoverable error in dgeev");
   else
     {
-      lambda.resize (n);
-      v.resize (n, n);
+      if (info > 0)
+	(*current_liboctave_error_handler) ("dgeev failed to converge");
+      else
+	{
+	  lambda.resize (n);
+	  v.resize (n, n);
 
-      for (int j = 0; j < n; j++)
-	{
-	  if (wi.elem (j) == 0.0)
+	  for (int j = 0; j < n; j++)
 	    {
-	      lambda.elem (j) = Complex (wr.elem (j));
-	      for (int i = 0; i < n; i++)
-		v.elem (i, j) = vr.elem (i, j);
-	    }
-	  else
-	    {
-	      if (j+1 >= n)
+	      if (wi.elem (j) == 0.0)
+		{
+		  lambda.elem (j) = Complex (wr.elem (j));
+		  for (int i = 0; i < n; i++)
+		    v.elem (i, j) = vr.elem (i, j);
+		}
+	      else
 		{
-		  (*current_liboctave_error_handler) ("EIG: internal error");
-		  return -1;
+		  if (j+1 >= n)
+		    {
+		      (*current_liboctave_error_handler)
+			("EIG: internal error");
+		      return -1;
+		    }
+
+		  for (int i = 0; i < n; i++)
+		    {
+		      lambda.elem(j) = Complex (wr.elem(j), wi.elem(j));
+		      lambda.elem(j+1) = Complex (wr.elem(j+1), wi.elem(j+1));
+		      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++;
 		}
+	    }
+	}
+    }
+
+  return info;
+}
+
+int
+EIG::symmetric_init (const Matrix& a)
+{
+  int n = a.rows ();
 
-	      for (int i = 0; i < n; i++)
-		{
-		  lambda.elem (j) = Complex (wr.elem (j), wi.elem (j));
-		  lambda.elem (j+1) = Complex (wr.elem (j+1), wi.elem (j+1));
-		  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++;
-	    }
+  if (n != a.cols ())
+    {
+      (*current_liboctave_error_handler) ("EIG requires square matrix");
+      return -1;
+    }
+
+  int info = 0;
+
+  Matrix atmp = a;
+  double *tmp_data = atmp.fortran_vec ();
+
+  Array<double> 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 ();
+
+  F77_XFCN (dsyev, DSYEV, ("V", "U", n, tmp_data, n, pwr, pwork,
+			   lwork, info, 1L, 1L));
+
+  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
+	{
+	  lambda.resize (n);
+
+	  for (int j = 0; j < n; j++)
+	    lambda.elem (j) = Complex (wr.elem (j));
+
+	  v = atmp;
 	}
     }
 
@@ -126,6 +195,9 @@
 int
 EIG::init (const ComplexMatrix& a)
 {
+  if (a.is_hermitian ())
+    return hermitian_init (a);
+
   int n = a.rows ();
 
   if (n != a.cols ())
@@ -134,33 +206,89 @@
       return -1;
     }
 
-  int info;
-
-  lambda.resize (n);
-  v.resize (n, n);
-  
-  Complex *pw = lambda.fortran_vec ();
-  Complex *pvr = v.fortran_vec ();
+  int info = 0;
 
   ComplexMatrix atmp = a;
   Complex *tmp_data = atmp.fortran_vec ();
 
+  ComplexColumnVector w (n);
+  Complex *pw = w.fortran_vec ();
+
+  ComplexMatrix vtmp (n, n);
+  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 ();
 
-  Array<double> rwork (4*n);
+  int lrwork = 2*n;
+  Array<double> rwork (lrwork);
   double *prwork = rwork.fortran_vec ();
 
   Complex *dummy = 0;
   int idummy = 1;
 
   F77_XFCN (zgeev, ZGEEV, ("N", "V", n, tmp_data, n, pw, dummy, idummy,
-			   pvr, n, pwork, lwork, prwork, info, 1L, 1L));
+			   pv, n, pwork, lwork, prwork, info, 1L, 1L));
+
+  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
+    {
+      lambda = w;
+      v = vtmp;
+    }
+
+  return info;
+}
+
+int
+EIG::hermitian_init (const ComplexMatrix& a)
+{
+  int n = a.rows ();
+
+  if (n != a.cols ())
+    {
+      (*current_liboctave_error_handler) ("EIG requires square matrix");
+      return -1;
+    }
+
+  int info = 0;
 
-  if (f77_exception_encountered)
-    (*current_liboctave_error_handler) ("unrecoverable error in zgeev");
+  ComplexMatrix atmp = a;
+  Complex *tmp_data = atmp.fortran_vec ();
+
+  ColumnVector w (n);
+  double *pw = w.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 lrwork = 3*n;
+  Array<double> rwork (lrwork);
+  double *prwork = rwork.fortran_vec ();
+
+  F77_XFCN (zheev, ZHEEV, ("V", "U", n, tmp_data, n, pw, pwork,
+			   lwork, prwork, info, 1L, 1L));
+
+  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
+    {
+      lambda = w;
+      v = atmp;
+    }
 
   return info;
 }
--- a/liboctave/EIG.h
+++ b/liboctave/EIG.h
@@ -80,6 +80,9 @@
 
   int init (const Matrix& a);
   int init (const ComplexMatrix& a);
+
+  int symmetric_init (const Matrix& a);
+  int hermitian_init (const ComplexMatrix& a);
 };
 
 #endif
--- a/readline/ChangeLog
+++ b/readline/ChangeLog
@@ -1,5 +1,8 @@
 Thu Mar 13 11:43:35 1997  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
+	* configure.in: Don't set special PIC options on SGI systems,
+	since PIC is apparently the default.
+
 	* Makefile.in (realclean, distclean): Don't remove configure or
 	config.h.in.
 
--- a/readline/configure.in
+++ b/readline/configure.in
@@ -49,6 +49,9 @@
     SHLEXT=sl
     SH_LDFLAGS="-shared -fPIC"
   ;;
+  *-sgi-*)
+    CPICFLAG=
+  ;;
 esac
 AC_MSG_RESULT([defining CPICFLAG to be $CPICFLAG])
 AC_MSG_RESULT([defining SHLEXT to be $SHLEXT])