--- a/liboctave/EIG.cc	Fri Mar 14 04:31:14 1997 +0000
+++ b/liboctave/EIG.cc	Fri Mar 14 08:25:58 1997 +0000
@@ -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;
 }