--- a/liboctave/ChangeLog	Fri Nov 21 15:03:03 2008 +0100
+++ b/liboctave/ChangeLog	Mon Nov 24 10:55:50 2008 +0100
@@ -1,3 +1,41 @@
+2008-11-21  Jarkko Kaleva  <d3roga@gmail.com>
+
+	* EIG.h (EIG::EIG (const Matrix& a, const Matrix& b, 
+	bool calc_eigenvectors = true)): New constructor.
+	(EIG::EIG (const Matrix& a, const Matrix& b, octave_idx_type& info, 
+	bool calc_eigenvectors = true)): New constructor.
+	(EIG::EIG (const ComplexMatrix& a, const ComplexMatrix& b, 
+	bool calc_eigenvectors = true)): New constructor.
+	(EIG::EIG (const ComplexMatrix& a, const ComplexMatrix& b, 
+	octave_idx_type& info, bool calc_eigenvectors = true)): New 
+	constructor.
+	* EIG.cc (EIG::init (const Matrix& a, const Matrix& b, 
+	bool calc_eigenvectors)): New function.
+	(EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, 
+	bool calc_eigenvectors)): New function.
+	(EIG::symmetric_init (const Matrix& a, const Matrix& b, 
+	bool calc_eigenvectors)): New function.
+	(EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, 
+	bool calc_eigenvectors)): New function.
+	* fEIG.h (fEIG::fEIG (const FloatMatrix& a, const FloatMatrix& b, 
+	bool calc_eigenvectors = true)): New constructor.
+	(fEIG::fEIG (const FloatMatrix& a, const FloatMatrix& b, 
+	octave_idx_type& info, bool calc_eigenvectors = true)): New 
+	constructor.
+	(fEIG::fEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b, 
+	bool calc_eigenvectors = true)): New constructor.	
+	(fEIG::fEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b, 
+	octave_idx_type& info, bool calc_eigenvectors = true)): New 
+	constructor.
+	(fEIG::init (const FloatMatrix& a, const FloatMatrix& b, 
+	bool calc_eigenvectors)): New function.
+	(fEIG::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, 
+	bool calc_eigenvectors)): New function.
+	(fEIG::symmetric_init (const FloatMatrix& a, const FloatMatrix& b, 
+	bool calc_eigenvectors)): New function.
+	(fEIG::hermitian_init (const FloatComplexMatrix& a, 
+	const FloatComplexMatrix& b, bool calc_eigenvectors)): New function.
+
 2008-11-19  Jaroslav Hajek  <highegg@gmail.com>
 	
 	* dMatrix.cc (Matrix::determinant),

--- a/liboctave/EIG.cc	Fri Nov 21 15:03:03 2008 +0100
+++ b/liboctave/EIG.cc	Mon Nov 24 10:55:50 2008 +0100
@@ -65,6 +65,68 @@
 			   Complex*, const octave_idx_type&, double*, octave_idx_type&
 			   F77_CHAR_ARG_LEN_DECL
 			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL, 
+			   const octave_idx_type&, 
+			   double*, const octave_idx_type&,
+			   octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, 
+			   const octave_idx_type&, 
+			   Complex*, const octave_idx_type&,
+			   octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dggev, DGGEV) (F77_CONST_CHAR_ARG_DECL, 
+			   F77_CONST_CHAR_ARG_DECL,
+			   const octave_idx_type&, 
+			   double*, const octave_idx_type&,
+			   double*, const octave_idx_type&,
+			   double*, double*, double *,
+			   double*, const octave_idx_type&, double*, const octave_idx_type&,
+			   double*, const octave_idx_type&, octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dsygv, DSYGV) (const octave_idx_type&,
+			   F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL,
+			   const octave_idx_type&, 
+			   double*, const octave_idx_type&,
+			   double*, const octave_idx_type&,
+			   double*, double*, const octave_idx_type&, octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zggev, ZGGEV) (F77_CONST_CHAR_ARG_DECL, 
+			   F77_CONST_CHAR_ARG_DECL,
+			   const octave_idx_type&, 
+			   Complex*, const octave_idx_type&,
+			   Complex*, const octave_idx_type&,
+			   Complex*, Complex*,
+			   Complex*, const octave_idx_type&, Complex*, const octave_idx_type&,
+			   Complex*, const octave_idx_type&, double*, octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zhegv, ZHEGV) (const octave_idx_type&,
+			   F77_CONST_CHAR_ARG_DECL, 
+			   F77_CONST_CHAR_ARG_DECL,
+			   const octave_idx_type&, 
+			   Complex*, const octave_idx_type&,
+			   Complex*, const octave_idx_type&,
+			   double*, Complex*, const octave_idx_type&, double*, 
+			   octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
 }
 
 octave_idx_type
@@ -391,6 +453,414 @@
   return info;
 }
 
+octave_idx_type
+EIG::init (const Matrix& a, const Matrix& b, bool calc_ev)
+{
+  if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
+    {
+      (*current_liboctave_error_handler)
+	("EIG: matrix contains Inf or NaN values");
+      return -1;
+    }
+
+  octave_idx_type n = a.rows ();
+  octave_idx_type nb = b.rows ();
+
+  if (n != a.cols () || nb != b.cols ())
+    {
+      (*current_liboctave_error_handler) ("EIG requires square matrix");
+      return -1;
+    }
+
+  if (n != nb)
+    {
+      (*current_liboctave_error_handler) ("EIG requires same size matrices");
+      return -1;
+    }
+
+  octave_idx_type info = 0;
+
+  Matrix tmp = b;
+  double *tmp_data = tmp.fortran_vec ();
+
+  F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+			     n, tmp_data, n, 
+			     info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+
+  if (a.is_symmetric () && b.is_symmetric () && info == 0)
+    return symmetric_init (a, b, calc_ev);
+
+  Matrix atmp = a;
+  double *atmp_data = atmp.fortran_vec ();
+
+  Matrix btmp = b;
+  double *btmp_data = btmp.fortran_vec ();
+
+  Array<double> ar (n);
+  double *par = ar.fortran_vec ();
+
+  Array<double> ai (n);
+  double *pai = ai.fortran_vec ();
+
+  Array<double> beta (n);
+  double *pbeta = beta.fortran_vec ();
+
+  volatile octave_idx_type nvr = calc_ev ? n : 0;
+  Matrix vr (nvr, nvr);
+  double *pvr = vr.fortran_vec ();
+
+  octave_idx_type lwork = -1;
+  double dummy_work;
+
+  double *dummy = 0;
+  octave_idx_type idummy = 1;
+
+  F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			   F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			   n, atmp_data, n, btmp_data, n, 
+			   par, pai, pbeta,
+			   dummy, idummy, pvr, n,
+			   &dummy_work, lwork, info
+			   F77_CHAR_ARG_LEN (1)
+			   F77_CHAR_ARG_LEN (1)));
+
+  if (info == 0)
+    {
+      lwork = static_cast<octave_idx_type> (dummy_work);
+      Array<double> work (lwork);
+      double *pwork = work.fortran_vec ();
+
+      F77_XFCN (dggev, DGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			       F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       n, atmp_data, n, btmp_data, n, 
+			       par, pai, pbeta,
+			       dummy, idummy, pvr, n,
+			       pwork, lwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+
+      if (info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in dggev");
+	  return info;
+	}
+
+      if (info > 0)
+	{
+	  (*current_liboctave_error_handler) ("dggev failed to converge");
+	  return info;
+	}
+
+      lambda.resize (n);
+      v.resize (nvr, nvr);
+
+      for (octave_idx_type j = 0; j < n; j++)
+	{
+	  if (ai.elem (j) == 0.0)
+	    {
+	      lambda.elem (j) = Complex (ar.elem (j) / beta.elem (j));
+	      for (octave_idx_type i = 0; i < nvr; i++)
+		v.elem (i, j) = vr.elem (i, j);
+	    }
+	  else
+	    {
+	      if (j+1 >= n)
+		{
+		  (*current_liboctave_error_handler) ("EIG: internal error");
+		  return -1;
+		}
+
+	      lambda.elem(j) = Complex (ar.elem(j) / beta.elem (j), 
+	                                ai.elem(j) / beta.elem (j));
+	      lambda.elem(j+1) = Complex (ar.elem(j+1) / beta.elem (j+1), 
+	                                  ai.elem(j+1) / beta.elem (j+1));
+
+	      for (octave_idx_type i = 0; i < nvr; i++)
+		{
+		  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++;
+	    }
+	}
+    }
+  else
+    (*current_liboctave_error_handler) ("dggev workspace query failed");
+
+  return info;
+}
+
+octave_idx_type 
+EIG::symmetric_init (const Matrix& a, const Matrix& b, bool calc_ev)
+{
+  octave_idx_type n = a.rows ();
+  octave_idx_type nb = b.rows ();
+
+  if (n != a.cols () || nb != b.cols ())
+    {
+      (*current_liboctave_error_handler) ("EIG requires square matrix");
+      return -1;
+    }
+
+  if (n != nb)
+    {
+      (*current_liboctave_error_handler) ("EIG requires same size matrices");
+      return -1;
+    }
+
+  octave_idx_type info = 0;
+
+  Matrix atmp = a;
+  double *atmp_data = atmp.fortran_vec ();
+
+  Matrix btmp = b;
+  double *btmp_data = btmp.fortran_vec ();
+
+  ColumnVector wr (n);
+  double *pwr = wr.fortran_vec ();
+
+  octave_idx_type lwork = -1;
+  double dummy_work;
+
+  F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			   F77_CONST_CHAR_ARG2 ("U", 1),
+			   n, atmp_data, n, 
+			   btmp_data, n, 
+			   pwr, &dummy_work, lwork, info
+			   F77_CHAR_ARG_LEN (1)
+			   F77_CHAR_ARG_LEN (1)));
+
+  if (info == 0)
+    {
+      lwork = static_cast<octave_idx_type> (dummy_work);
+      Array<double> work (lwork);
+      double *pwork = work.fortran_vec ();
+
+      F77_XFCN (dsygv, DSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       F77_CONST_CHAR_ARG2 ("U", 1),
+			       n, atmp_data, n, 
+			       btmp_data, n, 
+			       pwr, pwork, lwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+
+      if (info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in dsygv");
+	  return info;
+	}
+
+      if (info > 0)
+	{
+	  (*current_liboctave_error_handler) ("dsygv failed to converge");
+	  return info;
+	}
+
+      lambda = ComplexColumnVector (wr);
+      v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
+    }
+  else
+    (*current_liboctave_error_handler) ("dsygv workspace query failed");
+
+  return info;
+}
+
+octave_idx_type
+EIG::init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_ev)
+{
+  if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
+    {
+      (*current_liboctave_error_handler)
+	("EIG: matrix contains Inf or NaN values");
+      return -1;
+    }
+
+  octave_idx_type n = a.rows ();
+  octave_idx_type nb = b.rows ();
+
+  if (n != a.cols () || nb != b.cols())
+    {
+      (*current_liboctave_error_handler) ("EIG requires square matrix");
+      return -1;
+    }
+
+  if (n != nb)
+    {
+      (*current_liboctave_error_handler) ("EIG requires same size matrices");
+      return -1;
+    }
+
+  octave_idx_type info = 0;
+
+  ComplexMatrix tmp = b;
+  Complex*tmp_data = tmp.fortran_vec ();
+
+  F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+			     n, tmp_data, n, 
+			     info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+
+  if (a.is_hermitian () && b.is_hermitian () && info == 0)
+    return hermitian_init (a, calc_ev);
+
+  ComplexMatrix atmp = a;
+  Complex *atmp_data = atmp.fortran_vec ();
+
+  ComplexMatrix btmp = b;
+  Complex *btmp_data = btmp.fortran_vec ();
+
+  ComplexColumnVector alpha (n);
+  Complex *palpha = alpha.fortran_vec ();
+
+  ComplexColumnVector beta (n);
+  Complex *pbeta = beta.fortran_vec ();
+
+  octave_idx_type nvr = calc_ev ? n : 0;
+  ComplexMatrix vtmp (nvr, nvr);
+  Complex *pv = vtmp.fortran_vec ();
+
+  octave_idx_type lwork = -1;
+  Complex dummy_work;
+
+  octave_idx_type lrwork = 8*n;
+  Array<double> rwork (lrwork);
+  double *prwork = rwork.fortran_vec ();
+
+  Complex *dummy = 0;
+  octave_idx_type idummy = 1;
+
+  F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			   F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			   n, atmp_data, n, btmp_data, n, 
+			   palpha, pbeta, dummy, idummy,
+			   pv, n, &dummy_work, lwork, prwork, info
+			   F77_CHAR_ARG_LEN (1)
+			   F77_CHAR_ARG_LEN (1)));
+
+  if (info == 0)
+    {
+      lwork = static_cast<octave_idx_type> (dummy_work.real ());
+      Array<Complex> work (lwork);
+      Complex *pwork = work.fortran_vec ();
+
+      F77_XFCN (zggev, ZGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			       F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       n, atmp_data, n, btmp_data, n, 
+			       palpha, pbeta, dummy, idummy,
+			       pv, n, pwork, lwork, prwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+      
+      if (info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in zggev");
+	  return info;
+	}
+
+      if (info > 0)
+	{
+	  (*current_liboctave_error_handler) ("zggev failed to converge");
+	  return info;
+	}
+
+      lambda.resize (n);
+
+      for (octave_idx_type j = 0; j < n; j++)
+        lambda.elem (j) = alpha.elem (j) / beta.elem(j);
+
+      v = vtmp;
+    }
+  else
+    (*current_liboctave_error_handler) ("zggev workspace query failed");
+
+  return info;
+}
+
+octave_idx_type
+EIG::hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_ev)
+{
+  octave_idx_type n = a.rows ();
+  octave_idx_type nb = b.rows ();
+
+  if (n != a.cols () || nb != b.cols ())
+    {
+      (*current_liboctave_error_handler) ("EIG requires square matrix");
+      return -1;
+    }
+
+  if (n != nb)
+    {
+      (*current_liboctave_error_handler) ("EIG requires same size matrices");
+      return -1;
+    }
+
+  octave_idx_type info = 0;
+
+  ComplexMatrix atmp = a;
+  Complex *atmp_data = atmp.fortran_vec ();
+
+  ComplexMatrix btmp = b;
+  Complex *btmp_data = btmp.fortran_vec ();
+
+  ColumnVector wr (n);
+  double *pwr = wr.fortran_vec ();
+
+  octave_idx_type lwork = -1;
+  Complex dummy_work;
+
+  octave_idx_type lrwork = 3*n;
+  Array<double> rwork (lrwork);
+  double *prwork = rwork.fortran_vec ();
+
+  F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			   F77_CONST_CHAR_ARG2 ("U", 1),
+			   n, atmp_data, n, 
+			   btmp_data, n,
+			   pwr, &dummy_work, lwork,
+			   prwork, info
+			   F77_CHAR_ARG_LEN (1)
+			   F77_CHAR_ARG_LEN (1)));
+
+  if (info == 0)
+    {
+      lwork = static_cast<octave_idx_type> (dummy_work.real ());
+      Array<Complex> work (lwork);
+      Complex *pwork = work.fortran_vec ();
+
+      F77_XFCN (zhegv, ZHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       F77_CONST_CHAR_ARG2 ("U", 1),
+			       n, atmp_data, n, 
+			       btmp_data, n, 
+			       pwr, pwork, lwork, prwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+
+      if (info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in zhegv");
+	  return info;
+	}
+
+      if (info > 0)
+	{
+	  (*current_liboctave_error_handler) ("zhegv failed to converge");
+	  return info;
+	}
+
+      lambda = ComplexColumnVector (wr);
+      v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
+    }
+  else
+    (*current_liboctave_error_handler) ("zhegv workspace query failed");
+
+  return info;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***

--- a/liboctave/EIG.h	Fri Nov 21 15:03:03 2008 +0100
+++ b/liboctave/EIG.h	Mon Nov 24 10:55:50 2008 +0100
@@ -48,12 +48,24 @@
   EIG (const Matrix& a, octave_idx_type& info, bool calc_eigenvectors = true)
     { info = init (a, calc_eigenvectors); }
 
+  EIG (const Matrix& a, const Matrix& b, bool calc_eigenvectors = true)
+    { init (a, b, calc_eigenvectors); }
+
+  EIG (const Matrix& a, const Matrix& b, octave_idx_type& info, bool calc_eigenvectors = true)
+    { info = init (a, b, calc_eigenvectors); }
+
   EIG (const ComplexMatrix& a, bool calc_eigenvectors = true)
     { init (a, calc_eigenvectors); }
 
   EIG (const ComplexMatrix& a, octave_idx_type& info, bool calc_eigenvectors = true)
     { info = init (a, calc_eigenvectors); }
 
+  EIG (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_eigenvectors = true)
+    { init (a, b, calc_eigenvectors); }
+
+  EIG (const ComplexMatrix& a, const ComplexMatrix& b, octave_idx_type& info, bool calc_eigenvectors = true)
+    { info = init (a, b, calc_eigenvectors); }
+
   EIG (const EIG& a)
     : lambda (a.lambda), v (a.v) { }
 
@@ -81,10 +93,14 @@
   ComplexMatrix v;
 
   octave_idx_type init (const Matrix& a, bool calc_eigenvectors);
+  octave_idx_type init (const Matrix& a, const Matrix& b, bool calc_eigenvectors);
   octave_idx_type init (const ComplexMatrix& a, bool calc_eigenvectors);
+  octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_eigenvectors);
 
   octave_idx_type symmetric_init (const Matrix& a, bool calc_eigenvectors);
+  octave_idx_type symmetric_init (const Matrix& a, const Matrix& b, bool calc_eigenvectors);
   octave_idx_type hermitian_init (const ComplexMatrix& a, bool calc_eigenvectors);
+  octave_idx_type hermitian_init (const ComplexMatrix& a, const ComplexMatrix& b, bool calc_eigenvectors);
 };
 
 #endif

--- a/liboctave/fEIG.cc	Fri Nov 21 15:03:03 2008 +0100
+++ b/liboctave/fEIG.cc	Mon Nov 24 10:55:50 2008 +0100
@@ -65,6 +65,68 @@
 			   FloatComplex*, const octave_idx_type&, float*, octave_idx_type&
 			   F77_CHAR_ARG_LEN_DECL
 			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (spotrf, SPOTRF) (F77_CONST_CHAR_ARG_DECL, 
+			   const octave_idx_type&, 
+			   float*, const octave_idx_type&,
+			   octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (cpotrf, CPOTRF) (F77_CONST_CHAR_ARG_DECL, 
+			   const octave_idx_type&, 
+			   FloatComplex*, const octave_idx_type&,
+			   octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (sggev, SGGEV) (F77_CONST_CHAR_ARG_DECL, 
+			   F77_CONST_CHAR_ARG_DECL,
+			   const octave_idx_type&, 
+			   float*, const octave_idx_type&,
+			   float*, const octave_idx_type&,
+			   float*, float*, float*,
+			   float*, const octave_idx_type&, float*, const octave_idx_type&,
+			   float*, const octave_idx_type&, octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (ssygv, SSYGV) (const octave_idx_type&,
+			   F77_CONST_CHAR_ARG_DECL, F77_CONST_CHAR_ARG_DECL,
+			   const octave_idx_type&, 
+			   float*, const octave_idx_type&,
+			   float*, const octave_idx_type&,
+			   float*, float*, const octave_idx_type&, octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (cggev, CGGEV) (F77_CONST_CHAR_ARG_DECL, 
+			   F77_CONST_CHAR_ARG_DECL,
+			   const octave_idx_type&, 
+			   FloatComplex*, const octave_idx_type&,
+			   FloatComplex*, const octave_idx_type&,
+			   FloatComplex*, FloatComplex*,
+			   FloatComplex*, const octave_idx_type&, FloatComplex*, const octave_idx_type&,
+			   FloatComplex*, const octave_idx_type&, float*, octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (chegv, CHEGV) (const octave_idx_type&,
+			   F77_CONST_CHAR_ARG_DECL, 
+			   F77_CONST_CHAR_ARG_DECL,
+			   const octave_idx_type&, 
+			   FloatComplex*, const octave_idx_type&,
+			   FloatComplex*, const octave_idx_type&,
+			   float*, FloatComplex*, const octave_idx_type&, float*, 
+			   octave_idx_type&
+			   F77_CHAR_ARG_LEN_DECL
+			   F77_CHAR_ARG_LEN_DECL);
 }
 
 octave_idx_type
@@ -391,6 +453,414 @@
   return info;
 }
 
+octave_idx_type
+FloatEIG::init (const FloatMatrix& a, const FloatMatrix& b, bool calc_ev)
+{
+  if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
+    {
+      (*current_liboctave_error_handler)
+	("EIG: matrix contains Inf or NaN values");
+      return -1;
+    }
+
+  octave_idx_type n = a.rows ();
+  octave_idx_type nb = b.rows ();
+
+  if (n != a.cols () || nb != b.cols ())
+    {
+      (*current_liboctave_error_handler) ("EIG requires square matrix");
+      return -1;
+    }
+
+  if (n != nb)
+    {
+      (*current_liboctave_error_handler) ("EIG requires same size matrices");
+      return -1;
+    }
+
+  octave_idx_type info = 0;
+
+  FloatMatrix tmp = b;
+  float *tmp_data = tmp.fortran_vec ();
+
+  F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+			     n, tmp_data, n, 
+			     info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+
+  if (a.is_symmetric () && b.is_symmetric () && info == 0)
+    return symmetric_init (a, b, calc_ev);
+
+  FloatMatrix atmp = a;
+  float *atmp_data = atmp.fortran_vec ();
+
+  FloatMatrix btmp = b;
+  float *btmp_data = btmp.fortran_vec ();
+
+  Array<float> ar (n);
+  float *par = ar.fortran_vec ();
+
+  Array<float> ai (n);
+  float *pai = ai.fortran_vec ();
+
+  Array<float> beta (n);
+  float *pbeta = beta.fortran_vec ();
+
+  volatile octave_idx_type nvr = calc_ev ? n : 0;
+  FloatMatrix vr (nvr, nvr);
+  float *pvr = vr.fortran_vec ();
+
+  octave_idx_type lwork = -1;
+  float dummy_work;
+
+  float *dummy = 0;
+  octave_idx_type idummy = 1;
+
+  F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			   F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			   n, atmp_data, n, btmp_data, n, 
+			   par, pai, pbeta,
+			   dummy, idummy, pvr, n,
+			   &dummy_work, lwork, info
+			   F77_CHAR_ARG_LEN (1)
+			   F77_CHAR_ARG_LEN (1)));
+
+  if (info == 0)
+    {
+      lwork = static_cast<octave_idx_type> (dummy_work);
+      Array<float> work (lwork);
+      float *pwork = work.fortran_vec ();
+
+      F77_XFCN (sggev, SGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			       F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       n, atmp_data, n, btmp_data, n, 
+			       par, pai, pbeta,
+			       dummy, idummy, pvr, n,
+			       pwork, lwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+
+      if (info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in sggev");
+	  return info;
+	}
+
+      if (info > 0)
+	{
+	  (*current_liboctave_error_handler) ("sggev failed to converge");
+	  return info;
+	}
+
+      lambda.resize (n);
+      v.resize (nvr, nvr);
+
+      for (octave_idx_type j = 0; j < n; j++)
+	{
+	  if (ai.elem (j) == 0.0)
+	    {
+	      lambda.elem (j) = FloatComplex (ar.elem (j) / beta.elem (j));
+	      for (octave_idx_type i = 0; i < nvr; i++)
+		v.elem (i, j) = vr.elem (i, j);
+	    }
+	  else
+	    {
+	      if (j+1 >= n)
+		{
+		  (*current_liboctave_error_handler) ("EIG: internal error");
+		  return -1;
+		}
+
+	      lambda.elem(j) = FloatComplex (ar.elem(j) / beta.elem (j), 
+	                                     ai.elem(j) / beta.elem (j));
+	      lambda.elem(j+1) = FloatComplex (ar.elem(j+1) / beta.elem (j+1), 
+	                                       ai.elem(j+1) / beta.elem (j+1));
+
+	      for (octave_idx_type i = 0; i < nvr; i++)
+		{
+		  float real_part = vr.elem (i, j);
+		  float imag_part = vr.elem (i, j+1);
+		  v.elem (i, j) = FloatComplex (real_part, imag_part);
+		  v.elem (i, j+1) = FloatComplex (real_part, -imag_part);
+		}
+	      j++;
+	    }
+	}
+    }
+  else
+    (*current_liboctave_error_handler) ("sggev workspace query failed");
+
+  return info;
+}
+
+octave_idx_type 
+FloatEIG::symmetric_init (const FloatMatrix& a, const FloatMatrix& b, bool calc_ev)
+{
+  octave_idx_type n = a.rows ();
+  octave_idx_type nb = b.rows ();
+
+  if (n != a.cols () || nb != b.cols ())
+    {
+      (*current_liboctave_error_handler) ("EIG requires square matrix");
+      return -1;
+    }
+
+  if (n != nb)
+    {
+      (*current_liboctave_error_handler) ("EIG requires same size matrices");
+      return -1;
+    }
+
+  octave_idx_type info = 0;
+
+  FloatMatrix atmp = a;
+  float *atmp_data = atmp.fortran_vec ();
+
+  FloatMatrix btmp = b;
+  float *btmp_data = btmp.fortran_vec ();
+
+  FloatColumnVector wr (n);
+  float *pwr = wr.fortran_vec ();
+
+  octave_idx_type lwork = -1;
+  float dummy_work;
+
+  F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			   F77_CONST_CHAR_ARG2 ("U", 1),
+			   n, atmp_data, n, 
+			   btmp_data, n, 
+			   pwr, &dummy_work, lwork, info
+			   F77_CHAR_ARG_LEN (1)
+			   F77_CHAR_ARG_LEN (1)));
+
+  if (info == 0)
+    {
+      lwork = static_cast<octave_idx_type> (dummy_work);
+      Array<float> work (lwork);
+      float *pwork = work.fortran_vec ();
+
+      F77_XFCN (ssygv, SSYGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       F77_CONST_CHAR_ARG2 ("U", 1),
+			       n, atmp_data, n, 
+			       btmp_data, n, 
+			       pwr, pwork, lwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+
+      if (info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in ssygv");
+	  return info;
+	}
+
+      if (info > 0)
+	{
+	  (*current_liboctave_error_handler) ("ssygv failed to converge");
+	  return info;
+	}
+
+      lambda = FloatComplexColumnVector (wr);
+      v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
+    }
+  else
+    (*current_liboctave_error_handler) ("ssygv workspace query failed");
+
+  return info;
+}
+
+octave_idx_type
+FloatEIG::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool calc_ev)
+{
+  if (a.any_element_is_inf_or_nan () || b.any_element_is_inf_or_nan ())
+    {
+      (*current_liboctave_error_handler)
+	("EIG: matrix contains Inf or NaN values");
+      return -1;
+    }
+
+  octave_idx_type n = a.rows ();
+  octave_idx_type nb = b.rows ();
+
+  if (n != a.cols () || nb != b.cols())
+    {
+      (*current_liboctave_error_handler) ("EIG requires square matrix");
+      return -1;
+    }
+
+  if (n != nb)
+    {
+      (*current_liboctave_error_handler) ("EIG requires same size matrices");
+      return -1;
+    }
+
+  octave_idx_type info = 0;
+
+  FloatComplexMatrix tmp = b;
+  FloatComplex *tmp_data = tmp.fortran_vec ();
+
+  F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
+			     n, tmp_data, n, 
+			     info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+
+  if (a.is_hermitian () && b.is_hermitian () && info == 0)
+    return hermitian_init (a, calc_ev);
+
+  FloatComplexMatrix atmp = a;
+  FloatComplex *atmp_data = atmp.fortran_vec ();
+
+  FloatComplexMatrix btmp = b;
+  FloatComplex *btmp_data = btmp.fortran_vec ();
+
+  FloatComplexColumnVector alpha (n);
+  FloatComplex *palpha = alpha.fortran_vec ();
+
+  FloatComplexColumnVector beta (n);
+  FloatComplex *pbeta = beta.fortran_vec ();
+
+  octave_idx_type nvr = calc_ev ? n : 0;
+  FloatComplexMatrix vtmp (nvr, nvr);
+  FloatComplex *pv = vtmp.fortran_vec ();
+
+  octave_idx_type lwork = -1;
+  FloatComplex dummy_work;
+
+  octave_idx_type lrwork = 8*n;
+  Array<float> rwork (lrwork);
+  float *prwork = rwork.fortran_vec ();
+
+  FloatComplex *dummy = 0;
+  octave_idx_type idummy = 1;
+
+  F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			   F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			   n, atmp_data, n, btmp_data, n, 
+			   palpha, pbeta, dummy, idummy,
+			   pv, n, &dummy_work, lwork, prwork, info
+			   F77_CHAR_ARG_LEN (1)
+			   F77_CHAR_ARG_LEN (1)));
+
+  if (info == 0)
+    {
+      lwork = static_cast<octave_idx_type> (dummy_work.real ());
+      Array<FloatComplex> work (lwork);
+      FloatComplex *pwork = work.fortran_vec ();
+
+      F77_XFCN (cggev, CGGEV, (F77_CONST_CHAR_ARG2 ("N", 1),
+			       F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       n, atmp_data, n, btmp_data, n, 
+			       palpha, pbeta, dummy, idummy,
+			       pv, n, pwork, lwork, prwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+      
+      if (info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in cggev");
+	  return info;
+	}
+
+      if (info > 0)
+	{
+	  (*current_liboctave_error_handler) ("cggev failed to converge");
+	  return info;
+	}
+
+      lambda.resize (n);
+
+      for (octave_idx_type j = 0; j < n; j++)
+        lambda.elem (j) = alpha.elem (j) / beta.elem(j);
+
+      v = vtmp;
+    }
+  else
+    (*current_liboctave_error_handler) ("cggev workspace query failed");
+
+  return info;
+}
+
+octave_idx_type
+FloatEIG::hermitian_init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool calc_ev)
+{
+  octave_idx_type n = a.rows ();
+  octave_idx_type nb = b.rows ();
+
+  if (n != a.cols () || nb != b.cols ())
+    {
+      (*current_liboctave_error_handler) ("EIG requires square matrix");
+      return -1;
+    }
+
+  if (n != nb)
+    {
+      (*current_liboctave_error_handler) ("EIG requires same size matrices");
+      return -1;
+    }
+
+  octave_idx_type info = 0;
+
+  FloatComplexMatrix atmp = a;
+  FloatComplex *atmp_data = atmp.fortran_vec ();
+
+  FloatComplexMatrix btmp = b;
+  FloatComplex *btmp_data = btmp.fortran_vec ();
+
+  FloatColumnVector wr (n);
+  float *pwr = wr.fortran_vec ();
+
+  octave_idx_type lwork = -1;
+  FloatComplex dummy_work;
+
+  octave_idx_type lrwork = 3*n;
+  Array<float> rwork (lrwork);
+  float *prwork = rwork.fortran_vec ();
+
+  F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			   F77_CONST_CHAR_ARG2 ("U", 1),
+			   n, atmp_data, n, 
+			   btmp_data, n,
+			   pwr, &dummy_work, lwork,
+			   prwork, info
+			   F77_CHAR_ARG_LEN (1)
+			   F77_CHAR_ARG_LEN (1)));
+
+  if (info == 0)
+    {
+      lwork = static_cast<octave_idx_type> (dummy_work.real ());
+      Array<FloatComplex> work (lwork);
+      FloatComplex *pwork = work.fortran_vec ();
+
+      F77_XFCN (chegv, CHEGV, (1, F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
+			       F77_CONST_CHAR_ARG2 ("U", 1),
+			       n, atmp_data, n, 
+			       btmp_data, n, 
+			       pwr, pwork, lwork, prwork, info
+			       F77_CHAR_ARG_LEN (1)
+			       F77_CHAR_ARG_LEN (1)));
+
+      if (info < 0)
+	{
+	  (*current_liboctave_error_handler) ("unrecoverable error in zhegv");
+	  return info;
+	}
+
+      if (info > 0)
+	{
+	  (*current_liboctave_error_handler) ("zhegv failed to converge");
+	  return info;
+	}
+
+      lambda = FloatComplexColumnVector (wr);
+      v = calc_ev ? FloatComplexMatrix (atmp) : FloatComplexMatrix ();
+    }
+  else
+    (*current_liboctave_error_handler) ("zhegv workspace query failed");
+
+  return info;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***

--- a/liboctave/fEIG.h	Fri Nov 21 15:03:03 2008 +0100
+++ b/liboctave/fEIG.h	Mon Nov 24 10:55:50 2008 +0100
@@ -48,12 +48,24 @@
   FloatEIG (const FloatMatrix& a, octave_idx_type& info, bool calc_eigenvectors = true)
     { info = init (a, calc_eigenvectors); }
 
+  FloatEIG (const FloatMatrix& a, const FloatMatrix& b, bool calc_eigenvectors = true)
+    { init (a, b, calc_eigenvectors); }
+
+  FloatEIG (const FloatMatrix& a, const FloatMatrix& b, octave_idx_type& info, bool calc_eigenvectors = true)
+    { info = init (a, b, calc_eigenvectors); }
+
   FloatEIG (const FloatComplexMatrix& a, bool calc_eigenvectors = true)
     { init (a, calc_eigenvectors); }
 
   FloatEIG (const FloatComplexMatrix& a, octave_idx_type& info, bool calc_eigenvectors = true)
     { info = init (a, calc_eigenvectors); }
 
+  FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool calc_eigenvectors = true)
+    { init (a, b, calc_eigenvectors); }
+
+  FloatEIG (const FloatComplexMatrix& a, const FloatComplexMatrix& b, octave_idx_type& info, bool calc_eigenvectors = true)
+    { info = init (a, b, calc_eigenvectors); }
+
   FloatEIG (const FloatEIG& a)
     : lambda (a.lambda), v (a.v) { }
 
@@ -81,10 +93,14 @@
   FloatComplexMatrix v;
 
   octave_idx_type init (const FloatMatrix& a, bool calc_eigenvectors);
+  octave_idx_type init (const FloatMatrix& a, const FloatMatrix& b, bool calc_eigenvectors);
   octave_idx_type init (const FloatComplexMatrix& a, bool calc_eigenvectors);
+  octave_idx_type init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool calc_eigenvectors);
 
   octave_idx_type symmetric_init (const FloatMatrix& a, bool calc_eigenvectors);
+  octave_idx_type symmetric_init (const FloatMatrix& a, const FloatMatrix& b, bool calc_eigenvectors);
   octave_idx_type hermitian_init (const FloatComplexMatrix& a, bool calc_eigenvectors);
+  octave_idx_type hermitian_init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, bool calc_eigenvectors);
 };
 
 #endif

--- a/src/ChangeLog	Fri Nov 21 15:03:03 2008 +0100
+++ b/src/ChangeLog	Mon Nov 24 10:55:50 2008 +0100
@@ -1,3 +1,8 @@
+2008-11-21  Jarkko Kaleva  <d3roga@gmail.com>
+
+	* DLD-FUNCTIONS/eig.cc (Feig): Handle generalized eigenvalues and 
+	eigenvectors.
+
 2008-11-19  Jaroslav Hajek  <highegg@gmail.com>
 
 	* DLD_FUNCTIONS/det.cc: Include only DET.h. Retrieve & matrix type &

--- a/src/DLD-FUNCTIONS/eig.cc	Fri Nov 21 15:03:03 2008 +0100
+++ b/src/DLD-FUNCTIONS/eig.cc	Mon Nov 24 10:55:50 2008 +0100
@@ -37,7 +37,9 @@
 DEFUN_DLD (eig, args, nargout,
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {@var{lambda} =} eig (@var{a})\n\
+@deftypefnx {Loadable Function} {@var{lambda} =} eig (@var{a}, @var{b})\n\
 @deftypefnx {Loadable Function} {[@var{v}, @var{lambda}] =} eig (@var{a})\n\
+@deftypefnx {Loadable Function} {[@var{v}, @var{lambda}] =} eig (@var{a}, @var{b})\n\
 The eigenvalues (and eigenvectors) of a matrix are computed in a several\n\
 step process which begins with a Hessenberg decomposition, followed by a\n\
 Schur decomposition, from which the eigenvalues are apparent.  The\n\
@@ -51,55 +53,116 @@
 
   int nargin = args.length ();
 
-  if (nargin != 1 || nargout > 2)
+  if (nargin > 2 || nargin == 0 || nargout > 2)
     {
       print_usage ();
       return retval;
     }
 
-  octave_value arg = args(0);
+  octave_value arg_a, arg_b;
+
+  octave_idx_type nr_a = 0, nr_b = 0;
+  octave_idx_type nc_a = 0, nc_b = 0;
 
-  octave_idx_type nr = arg.rows ();
-  octave_idx_type nc = arg.columns ();
+  arg_a = args(0);
+  nr_a = arg_a.rows ();
+  nc_a = arg_a.columns ();
 
-  int arg_is_empty = empty_arg ("eig", nr, nc);
+  int arg_is_empty = empty_arg ("eig", nr_a, nc_a);
   if (arg_is_empty < 0)
     return retval;
   else if (arg_is_empty > 0)
     return octave_value_list (2, Matrix ());
 
-  if (nr != nc)
+  if (!(arg_a.is_single_type () || arg_a.is_double_type ()))
+    {
+      gripe_wrong_type_arg ("eig", arg_a);
+      return retval;
+    }
+
+  if (nargin == 2)
+    {
+      arg_b = args(1);
+      nr_b = arg_b.rows ();
+      nc_b = arg_b.columns ();
+
+      arg_is_empty = empty_arg ("eig", nr_b, nc_b);
+      if (arg_is_empty < 0)
+        return retval;
+      else if (arg_is_empty > 0)
+        return octave_value_list (2, Matrix ());
+
+      if (!(arg_b.is_single_type() || arg_b.is_double_type ()))
+	{
+	  gripe_wrong_type_arg ("eig", arg_b);
+	  return retval;
+	}
+    }
+
+  if (nr_a != nc_a)
     {
       gripe_square_matrix_required ("eig");
       return retval;
     }
 
-  Matrix tmp;
-  ComplexMatrix ctmp;
-  FloatMatrix ftmp;
-  FloatComplexMatrix fctmp;
+  if (nargin == 2 && nr_b != nc_b)
+    {
+      gripe_square_matrix_required ("eig");
+      return retval;
+    }
 
-  if (arg.is_single_type ())
+  Matrix tmp_a, tmp_b;
+  ComplexMatrix ctmp_a, ctmp_b;
+  FloatMatrix ftmp_a, ftmp_b;
+  FloatComplexMatrix fctmp_a, fctmp_b;
+
+  if (arg_a.is_single_type ())
     {
       FloatEIG result;
 
-      if (arg.is_real_type ())
+      if (nargin == 1)
 	{
-	  ftmp = arg.float_matrix_value ();
+	  if (arg_a.is_real_type ())
+	    {
+	      ftmp_a = arg_a.float_matrix_value ();
 
-	  if (error_state)
-	    return retval;
+	      if (error_state)
+	        return retval;
+	      else
+	        result = FloatEIG (ftmp_a, nargout > 1);
+	    }
 	  else
-	    result = FloatEIG (ftmp, nargout > 1);
+	    {
+	      fctmp_a = arg_a.float_complex_matrix_value ();
+
+	      if (error_state)
+	        return retval;
+	      else
+	        result = FloatEIG (fctmp_a, nargout > 1);
+	    }
 	}
-      else if (arg.is_complex_type ())
+      else if (nargin == 2)
 	{
-	  fctmp = arg.float_complex_matrix_value ();
+	  if (arg_a.is_real_type () && arg_b.is_real_type ())
+	    {
+	      ftmp_a = arg_a.float_matrix_value ();
+	      ftmp_b = arg_b.float_matrix_value ();
 
-	  if (error_state)
-	    return retval;
+	      if (error_state)
+	        return retval;
+	      else
+	        result = FloatEIG (ftmp_a, ftmp_b, nargout > 1);
+	    }
 	  else
-	    result = FloatEIG (fctmp, nargout > 1);
+	    {
+	      fctmp_a = arg_a.float_complex_matrix_value ();
+	      fctmp_b = arg_b.float_complex_matrix_value ();
+
+	      if (error_state)
+	        return retval;
+	      else
+	        result = FloatEIG (fctmp_a, fctmp_b, nargout > 1);
+	    }
 	}
 
       if (! error_state)
@@ -123,28 +186,49 @@
     {
       EIG result;
 
-      if (arg.is_real_type ())
-	{
-	  tmp = arg.matrix_value ();
-
-	  if (error_state)
-	    return retval;
-	  else
-	    result = EIG (tmp, nargout > 1);
-	}
-      else if (arg.is_complex_type ())
+      if (nargin == 1)
 	{
-	  ctmp = arg.complex_matrix_value ();
+	  if (arg_a.is_real_type ())
+	    {
+	      tmp_a = arg_a.matrix_value ();
 
-	  if (error_state)
-	    return retval;
+	      if (error_state)
+	        return retval;
+	      else
+	        result = EIG (tmp_a, nargout > 1);
+	    }
 	  else
-	    result = EIG (ctmp, nargout > 1);
+	    {
+	      ctmp_a = arg_a.complex_matrix_value ();
+
+	      if (error_state)
+	        return retval;
+	      else
+	        result = EIG (ctmp_a, nargout > 1);
+	    }
 	}
-      else
+      else if (nargin == 2)
 	{
-	  gripe_wrong_type_arg ("eig", tmp);
-	  return retval;
+	  if (arg_a.is_real_type () && arg_b.is_real_type ())
+	    {
+	      tmp_a = arg_a.matrix_value ();
+	      tmp_b = arg_b.matrix_value ();
+
+	      if (error_state)
+	        return retval;
+	      else
+	        result = EIG (tmp_a, tmp_b, nargout > 1);
+	    }
+	  else 
+	    {
+	      ctmp_a = arg_a.complex_matrix_value ();
+	      ctmp_b = arg_b.complex_matrix_value ();
+
+	      if (error_state)
+	        return retval;
+	      else
+	        result = EIG (ctmp_a, ctmp_b, nargout > 1);
+	    }
 	}
 
       if (! error_state)
@@ -186,9 +270,67 @@
 %! assert(d, single([-1, 0; 0, 3]), sqrt (eps('single')));
 %! assert(v, [-x, x; x, x], sqrt (eps('single')));
 
+%!test
+%! A = [1, 2; -1, 1]; B = [3, 3; 1, 2];
+%! [v, d] = eig (A, B);
+%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
+%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
+
+%!test
+%! A = single([1, 2; -1, 1]); B = single([3, 3; 1, 2]);
+%! [v, d] = eig (A, B);
+%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps('single')));
+%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps('single')));
+
+%!test
+%! A = [1, 2; 2, 1]; B = [3, -2; -2, 3];
+%! [v, d] = eig (A, B);
+%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
+%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
+
+%!test
+%! A = single([1, 2; 2, 1]); B = single([3, -2; -2, 3]);
+%! [v, d] = eig (A, B);
+%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps('single')));
+%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps('single')));
+
+%!test
+%! A = [1+3i, 2+i; 2-i, 1+3i]; B = [5+9i, 2+i; 2-i, 5+9i];
+%! [v, d] = eig (A, B);
+%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
+%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
+
+%!test
+%! A = single([1+3i, 2+i; 2-i, 1+3i]); B = single([5+9i, 2+i; 2-i, 5+9i]);
+%! [v, d] = eig (A, B);
+%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps('single')));
+%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps('single')));
+
+%!test
+%! A = [1+3i, 2+3i; 3-8i, 8+3i]; B = [8+i, 3+i; 4-9i, 3+i];
+%! [v, d] = eig (A, B);
+%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
+%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
+
+%!test
+%! A = single([1+3i, 2+3i; 3-8i, 8+3i]); B = single([8+i, 3+i; 4-9i, 3+i]);
+%! [v, d] = eig (A, B);
+%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps('single')));
+%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps('single')));
+
+%!test
+%! A = [1, 2; 3, 8]; B = [8, 3; 4, 3];
+%! [v, d] = eig (A, B);
+%! assert(A * v(:, 1), d(1, 1) * B * v(:, 1), sqrt (eps));
+%! assert(A * v(:, 2), d(2, 2) * B * v(:, 2), sqrt (eps));
+
 %!error <Invalid call to eig.*> eig ();
-%!error <Invalid call to eig.*> eig ([1, 2; 3, 4], 2);
+%!error <Invalid call to eig.*> eig ([1, 2; 3, 4], [4, 3; 2, 1], 1);
+%!error eig ([1, 2; 3, 4], 2);
 %!error eig ([1, 2; 3, 4; 5, 6]);
+%!error eig ("abcd");
+%!error eig ([1 2 ; 2 3], "abcd");
+%!error eig (false, [1 2 ; 2 3]);
 
  */