Generalized eigenvalues

[Jarkko Kaleva]

Hi,

I decided to take a shot at implementing a feature for solving
generalized eigenvalue problems.

For EIG and fEIG classes the implementation was quite straightforward.
Added the needed lapack calls, otherwise the functions are very close
to the basic eigenvalue functions.

DLD-FUNCTIONS/eig.cc needed a bit more work. But I think that there
wasn't too much increase in complexity. The original structure is
still intact, except for the matrix type check, which is now done
before calculating the eigenvalues.

I added test cases for following input cases (single and double value):
  - Nonsymmetric real, nonsymmetric real
  - Symmetric real, symmetric and positive definite real
  - Hermitian complex, hermitian positive definite complex
  - Non-hermitian complex, non-hermitian complex
These tests should go through all the additions to the EIG and fEIG
classes. (Maybe these should be located in EIG.cc and fEIG.cc)

I also added a couple tests to check the input of invalid matrix types.

-- Jarkko Kaleva

[generalized-eigenvalues.patch]

---

diff -ru a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog 2008-08-21 17:19:13.000000000 +0300
+++ b/liboctave/ChangeLog 2008-08-22 18:33:56.000000000 +0300
@@ -1,3 +1,42 @@
[phone]  Jarkko Kaleva  <d3roga@...>
+
+ * 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-08-19  David Bateman  <dbateman@...>
 
  * oct-inttypes.h (template <class T1, class T2> inline T2
diff -ru a/liboctave/EIG.cc b/liboctave/EIG.cc
--- a/liboctave/EIG.cc 2008-08-21 17:19:13.000000000 +0300
+++ b/liboctave/EIG.cc 2008-08-22 18:33:56.000000000 +0300
@@ -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++ ***
diff -ru a/liboctave/EIG.h b/liboctave/EIG.h
--- a/liboctave/EIG.h 2008-08-21 17:19:13.000000000 +0300
+++ b/liboctave/EIG.h 2008-08-22 18:33:56.000000000 +0300
@@ -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
diff -ru a/liboctave/fEIG.cc b/liboctave/fEIG.cc
--- a/liboctave/fEIG.cc 2008-08-21 17:19:13.000000000 +0300
+++ b/liboctave/fEIG.cc 2008-08-22 18:33:56.000000000 +0300
@@ -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++ ***
diff -ru a/liboctave/fEIG.h b/liboctave/fEIG.h
--- a/liboctave/fEIG.h 2008-08-21 17:19:13.000000000 +0300
+++ b/liboctave/fEIG.h 2008-08-22 18:33:56.000000000 +0300
@@ -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
diff -ru a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog 2008-08-22 16:37:25.000000000 +0300
+++ b/src/ChangeLog 2008-08-22 18:33:56.000000000 +0300
@@ -1,3 +1,8 @@
[phone]  Jarkko Kaleva  <d3roga@...>
+
+ * DLD-FUNCTIONS/eig.cc (Feig): Handle generalized eigenvalues and
+ eigenvectors.
+
 2008-08-21  Thomas Treichl  <Thomas.Treichl@...>
 
  * mappers.cc: Increase test script tolerance.
diff -ru a/src/DLD-FUNCTIONS/eig.cc b/src/DLD-FUNCTIONS/eig.cc
--- a/src/DLD-FUNCTIONS/eig.cc 2008-08-21 17:19:15.000000000 +0300
+++ b/src/DLD-FUNCTIONS/eig.cc 2008-08-22 18:33:56.000000000 +0300
@@ -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 = arg.rows ();
-  octave_idx_type nc = arg.columns ();
+  octave_idx_type nr_a = 0, nr_b = 0;
+  octave_idx_type nc_a = 0, nc_b = 0;
 
-  int arg_is_empty = empty_arg ("eig", nr, nc);
+  arg_a = args(0);
+  nr_a = arg_a.rows ();
+  nc_a = arg_a.columns ();
+
+  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;
+    }
+
+  if (nargin == 2 && nr_b != nc_b)
     {
       gripe_square_matrix_required ("eig");
       return retval;
     }
 
-  Matrix tmp;
-  ComplexMatrix ctmp;
-  FloatMatrix ftmp;
-  FloatComplexMatrix fctmp;
+  Matrix tmp_a, tmp_b;
+  ComplexMatrix ctmp_a, ctmp_b;
+  FloatMatrix ftmp_a, ftmp_b;
+  FloatComplexMatrix fctmp_a, fctmp_b;
 
-  if (arg.is_single_type ())
+  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 ())
+      if (nargin == 1)
  {
-  tmp = arg.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 (tmp, nargout > 1);
- }
-      else if (arg.is_complex_type ())
- {
-  ctmp = arg.complex_matrix_value ();
+    {
+      ctmp_a = arg_a.complex_matrix_value ();
 
-  if (error_state)
-    return retval;
-  else
-    result = EIG (ctmp, nargout > 1);
+      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]);
 
  */
 

---

[Khaled]

Dear Jarkko Kaleva,

   I am really in need for the generalized eigenvector on the form eig(A,B) and I liked your patch a lot. However, I am trying to apply this patch on liboctave/EIG.cc, EIG.h and  src/DLD-FUNCTIONS/eig.cc. The later does not work, and I have manually figured out that the file is different from your in the patch. I always get this error:

patching file src/DLD-FUNCTIONS/eig.cc
Hunk #1 succeeded at 36 (offset -1 lines).
Hunk #2 FAILED at 52.
Hunk #3 FAILED at 185.
Hunk #4 FAILED at 269.
3 out of 4 hunks FAILED -- saving rejects to file src/DLD-FUNCTIONS/eig.cc.rej


I am using octave 3.0.3, can you please help?

 

[Jarkko Kaleva-2]

Hi,

After that post I have revised the patch once, and it is now added to
the development sources.
So the easiest way to get the feature might be by compiling octave
from the development sources.

On Sat, Jan 3, 2009 at 11:40 PM, Khaled <kshawky@...> wrote:

>
> Dear Jarkko Kaleva,
>
>   I am really in need for the generalized eigenvector on the form eig(A,B)
> and I liked your patch a lot. However, I am trying to apply this patch on
> liboctave/EIG.cc, EIG.h and  src/DLD-FUNCTIONS/eig.cc. The later does not
> work, and I have manually figured out that the file is different from your
> in the patch. I always get this error:
>
> patching file src/DLD-FUNCTIONS/eig.cc
> Hunk #1 succeeded at 36 (offset -1 lines).
> Hunk #2 FAILED at 52.
> Hunk #3 FAILED at 185.
> Hunk #4 FAILED at 269.
> 3 out of 4 hunks FAILED -- saving rejects to file
> src/DLD-FUNCTIONS/eig.cc.rej
>
>
> I am using octave 3.0.3, can you please help?
>
>
>
> --
> View this message in context: http://www.nabble.com/Generalized-eigenvalues-tp19109913p21269754.html
> Sent from the Octave - Maintainers mailing list archive at Nabble.com.
>
>

--
Jarkko Kaleva