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 @@
+2008-08-22 Jarkko Kaleva
+
+ * 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
* oct-inttypes.h (template 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 ar (n);
+ double *par = ar.fortran_vec ();
+
+ Array ai (n);
+ double *pai = ai.fortran_vec ();
+
+ Array 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 (dummy_work);
+ Array 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 (dummy_work);
+ Array 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 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 (dummy_work.real ());
+ Array 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 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 (dummy_work.real ());
+ Array 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 ar (n);
+ float *par = ar.fortran_vec ();
+
+ Array ai (n);
+ float *pai = ai.fortran_vec ();
+
+ Array 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 (dummy_work);
+ Array 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 (dummy_work);
+ Array 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 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 (dummy_work.real ());
+ Array 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 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 (dummy_work.real ());
+ Array 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 @@
+2008-08-22 Jarkko Kaleva
+
+ * DLD-FUNCTIONS/eig.cc (Feig): Handle generalized eigenvalues and
+ eigenvectors.
+
2008-08-21 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} address@hidden =} eig (@var{a})\n\
address@hidden {Loadable Function} address@hidden =} eig (@var{a}, @var{b})\n\
@deftypefnx {Loadable Function} address@hidden, @var{lambda}] =} eig (@var{a})\n\
address@hidden {Loadable Function} address@hidden, @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 eig ();
-%!error eig ([1, 2; 3, 4], 2);
+%!error 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]);
*/